Ind. Eng. Chem. Res. 1995,34, 1219-1227
1219
PROCESS DESIGN AND CONTROL On-Line Estimation of Reaction Rates in Semicontinuous Reactors Guillermo E. Eliqabe3 Eser Ozdeger, and Christos Georgakis" Chemical Process Modeling and Control Research Center and Department of Chemical Engineering, Lehigh University, Bethlehem, Pennsylvania 18015
Cajetan Cordeiro Air Products and Chemicals Inc., Allentown, Pennsylvania 18195
State estimation techniques are used to on-line estimate reaction rates in well-mixed semicontinuous reactors. The proposed scheme assumes the availability of on-line measurements of the total mass of the reactants in the reactor, and copes with the lack of knowledge of the reaction rate dependence with concentration. Concepts from system theory are used to reduce the degrees of freedom that the pure Kalman filtering design produces. By the use of three examples of semicontinuous emulsion copolymerization, it is shown that accurate estimates of the reaction rates are possible in the presence of inaccurate measurements of the incoming flow rates and much earlier than the establishment of the pseudo-steady-state condition. The proposed technique is applied t o the experimental data collected during the semibatch polymerization of the following ternary systems: (a) acrylonitrile-ethyl acrylate-N-methylacrylamide and (b) butyl acrylate-styrene-acrylic acid. 1. Introduction
Very often one needs to know the current values of reaction rates in order to control chemical reactors effectively or in order to avoid unsafe operating conditions. When these reaction rates vary with time, as is the case in discontinuous reactors, the need to know how they evolve becomes more important. A typical example of this is the polymerization of two or more monomers in a discontinuous reactor. There, the composition of the copolymer formed during the reaction, in terms of one of the monomers, will be determined at any instant by the ratio of the rate at which that monomer reacts, to the total rate of polymerization (Ray and Gall, 1969; Tirrell and Gromley, 1981). This quantity, called instantaneous copolymer composition, plays an important role in determining many of the final properties of the produced copolymer. If a perfectly accurate model of the process were available, the values of these reaction rates could be calculated on-line and in parallel with the process operation. However, such an accurate model is rarely available, mostly because, in many reacting systems of practical interest, the reaction rates are complex functions of the concentrations of the reacting species, and the kinetic parameters involved are several. For example in heterogeneous systems, such as emulsion polymerizations, rate expressions do not depend on overall concentrations but are functions of the concentrations in a particular phase (Min and Ray, 1974).This
* To whom correspondence should be addressed: Chemical Process Modeling and Control Research Center, 111 Research Drive, Lehigh University, Bethlehem, PA 18015-4732. Phone: (610) 758-5432; Fax: (610) 758-5297; E-mail: cgOO@ Lehigh.edu. Present Address: Departamento de Ingenierfa Quimica, Universidad Nacional de Mar del Plata, J.B. Just0 4302, (7600) Mar Del Plata, Argentina. +
0888-588519512634-1219$09.00/0
then further necessitates the understanding of the thermodynamic equilibrium characteristics between the different phases and increases the difficulty of model development. In this paper we will show that even when no detailed model of the process is available it is still possible to obtain on-line estimates of reaction rates from on-line measurements of the total amounts of participating species introduced inside the reactor. This is accomplished using a very simple model of the process; this model is basically derived from the definition of rate of reaction and is independent of the type of reaction that is being carried out. Because the resulting equations for each component of this simple model are independent from each other, each reaction rate can be independently calculated from the rest. If a component cannot be measured, reaction rate estimates can still be calculated for those that can be measured. The estimation scheme will be based on Kalman filtering techniques (Stengel, 1986; Lewis, 1986). In the past other authors addressed the problem of estimating reaction rates using Kalman filter based algorithms (Schuler and Schmidt, 1992). These authors analyzed the use of calorimetric balances in the estimation of reaction rates, conversion, and rate of production, and pointed out the importance of estimating these variables to control runaways in batch reactors and to avoid overfeeding and runaways in semibatch reactors. In the following sections, rigorous results will be presented for two types of estimators, the continuous and the discrete one. The more intuitive but less realistic continuous estimator, assuming the availability of continuous measurements, will be used in the main part of the derivation to specify all the parameters needed for implementation. Because concentration measurements are necessarily discrete, the results obtained for the continuous case will be used to design the discrete type of the estimator with the same
0 1995 American Chemical Society
1220 Ind. Eng. Chem. Res., Vol. 34,No. 4,1995
characteristics already specified for the continuous design. Both the continuous and discrete algorithms will be tested with simulated data obtained from a very detailed model of the copolymerization of vinyl acetate and butyl acrylate (Dimitratos et al., 1991). Then, experimental data obtained during the copolymerization of acrylonitrile and ethyl acrylate as well as data from the copolymerization of styrene and butyl acrylate will be used to test the discrete version of the estimator.
with k1 = p d r and k2 = p1dr as the estimator gains. pI1andp12 are respectively the (1,l)and (1,2) elements of the (2 x 2) state covariance matrix whose timevarying elements are described by the following Riccati differential equations:
--
(7)
dt
dP22 2. Theory
The following n equations describe the evolution in time ( t )of the total mass (Mi) of the n different species participating in a semicontinuous reaction, dMildt = Fi - Ri for z = 1,..., ni
(1)
The chemical species are added t o the reactor a t a mass flow rate Fi and consumed by chemical reaction at a mass reaction rate Ri. Equations 1 allow the treatment of the particular case in which Fi = 0 for all i, which corresponds to a batch reaction. The fair assumption that the flow rates are known without error is made and the common case in which one has no knowledge of the kinetic expressions is considered. Under this last assumption eqs 1 become uncoupled and then only one equation a t a time needs to be considered to develop the respective reaction rate estimator. Although this simple equation is absolutely correct by definition, it is also of little use if nothing is said about Ri. In order to cast the model equations in terms of the stochastic estimation theory, and because no knowledge of the reaction rate is available, the simple random walk model is assumed for the Ri variables. Thus the final reactor model for a generic species is dM1dt = F - R
+ w1
r + 92
dt dp12 - -p22 --
M,=M+v
(4)
where the subscript i was dropped for simplicity. Equation 4 expresses how the measurement (M,) of the total mass of a generic component in the reactor is related to the system states. The standard assumptions that w1, w2, and u are white zero mean Gaussian random processes with spectral densities 41, q 2 , and r are made. It is also assumed that wl, w2, and u are uncorrelated. Under the previous assumptions eq 3 is the random walk model for the reaction rate as mentioned above. As can be seen, the input w1 was added to the first equation. As we will see later, this gives greater flexibility when specifying the desired performance of the estimator. Notice that because eqs 2-4 are linear, a standard linear Kalman filter to estimate M and R can be easily developed. 3. The Continuous Estimator
Using the standard Kalman filter algorithm (Schuler and Schmidt, 1992) the estimator equations can be written as &ldt
=F -R
+ kl(Mm -
&ldt = k2(Mm- hfj
(5)
(6)
r
(9)
Because the covariance matrix is symmetrical, just three equations describe its evolution. These equations depend on three parameters: 41, q2, and r. However from eqs 7-9 it is clear that the evolution of k l and k2 depend only on the ratios qllr and qdr. Thus two parameters need t o be chosen in order to completely determine the estimator. In many practical applications these parameters are used as tuning parameters and are selected by trial and error. In what follows we will show that, for the estimation problem that is being addressed, there is no need to choose these parameters by trial and error. For the case of a fixed gain Kalman filter, the use of linear analysis and stochastic systems theory tools will allow us to determine all the “adjustable” parameters based on our knowledge of the system and the type of behavior that we are expecting for the estimator. The steady state gain for the fxed-gain Kalman filter is calculated using the only stable steady state of eqs 7-9. For that steady state
(2)
(3)
- PllPl2 -
dt
k aRldt = w2
Pl22
=-e)
(10)
112
2
After setting the gains to constant values, the equations of the estimator (eqs 5 and 6)become linear and can be analyzed using Lqplace techniques. Thus, the estimated reaction rate R(s) can be expressed in the Laplace domain as
The roots of the denominator of eq 12 will determine the behavior of the estimator. As said before, these roots depend only on gllr and gdr. As qllr varies from -2(qdr)ll2 to 2(q2/rP2,the two roots take values from ~t$qdr)~/~ to (q2/r)1/4,describing a semicircle with radius (q2/r)u4on the s-plane [see Figure 11. If qllr is further increased, one of the roots tends to --m and the other one tends to 0. The estimator is stable provided qllr -2(qd~-)l/~ and qdr > 0. If q1lr is set to 0, the roots of the denominator polynomial are a complex conjugate pair located on lines forming angles of f 4 5 ” with the real negative axis of the s-plane. The damping ratio in this case is 4212, which corresponds t o the “optimal” damping ratio that provides minimum settling time for a second order system. Although a value of qllr different from 0 can be selected, in most cases 0 is an appropriate and realistic choice. This is primarily
Ind. Eng. Chem. Res., Vol. 34, No. 4, 1995 1221 systems theory (Papoulis 1984, p 3131,
where E represents expected value and y ( t )is the output of a linear system [with impulse response h(t)l driven by white noise of spectral density r. It is possible t o obtain a relationship between the parameter to be determined qdr and the tolerable amount 0: noise in the estimation. Therefore, by setting y ( t ) R(z)(t)and calculating the impulse response h(t) of the linear system defined by the lower block in Figure 2, it can be %hewn that the variance of the noisy part of R(t),i.e., R(2)(t), is given by
p = (&14)(qdr)~’~r Figure 1. Locus of the roots of the denominator polynomial of the transfer function representation of the estimator as qllr is varied from -2(q2/r)1/2 to -.
Figure 2. Equivalent transfer function representation of the estimated rate of reaction.
justifiable because eq 2 is, in principle, the definition of the reaction rate R provided that all incoming and outgoing flow rates have been accounted in the values of F. Now, it only remains to set the value of qdr that will determine the size of the pair of eigenvalues that will govern the behavior of the estimator. If the Laplace representation of the estimator (eq 12) is rearranged as shown in Figure 2, the estimated reaction rate is now represented*as the sum of two signals. On one hand, the signal R%) will contain the desired information smoothed by a second order filter; the smoothing effect should be as little as possible in order to preserve the identitx of the original signal. On the other hand, the signal R(2)(s) will not have any useful information and should be kept as small as possible. Clearly these are two conflicting objectives. As the smoothing effect is reduced, the noise is amplified. However it is always possible to specify the worst admissible condition for one of the disturbing effects (either smoothing or noise amplification) and design the estimator to meet that specification. If as a result the effect not considered in the design produces an inadmissible behavior, the requirement imposed a t the design stage should be relaxed until an “optimum” estimator is obtained in which a compromise between speed (or tracking capability) and level of noise is obtained. For simplicity, the level of noise that can be accepted will be used as the determining factor in the selection of the remaining parameter qdr that still needs to be chosen. The “colored” noise R(2)(s)that will be constrained to a certain range of admissible values is the result of white noise filtered through a second order filter as shown in Figure 2. Using the following result from stochastic
(14)
Thus a value for qdr can be calculated using the spectral density r of the measurement noise u (normally an estimate of this value is available) and the “desired variance p for the output noise R(2). The smaller the variance of this “colored” output noise is chosen, the slower the estimator response becomes. As mentioned above, there will be normally a trade off between filter speed and amount of noise in the estimation. The guidelines to design the state estimator given in the previous paragraphs are based on the availability of continuous measurements. However, that is not normally the case in the applications we are discussing here. Therefore a discrete version of the estimator is desirable. In order to take advantage of the guidelines developed for the more familiar continuous estimator and make use of the values of qllr and qdr obtained under those guidelines, the continuous model of eqs 2-4 will be discretized and covariances Q d l l , qdlz = qdzl (in the discrete version the sampling process correlates the originally uncorrelated noises w1 and WZ),q d Z 2 , and rd, equivalent to the spectral densities q1,qz, and r, will be found. On the basis of these values an optimum discrete Kalman filter for the discretized process will be designed. 4. The Discrete Estimator
A discrete version of the continuous model given by eqs 2-4 can be readily developed because this model is linear. If F is piecewise continuous and the sampling time is constant, the sampled data equations equivalent to the continuous model are
M ( k + l ) = M ( k ) - R(k)At + F(k)At + w , ( k )
+ = M ( k ) + v(k)
R ( k + l ) = R(k) ulz(k)
M&)
(15) (16) (17)
Now instead of white zero mean Gaussian random processes, the variables W I ,wg, and u are white zero mean Gaussian random sequences whose covariance matrices are numerically different from the corresponding spectral density matrices of the “continuous”noises. Thus, in order to have a rigorous sample data representation of the continuous model, the covariance matrices must be calculated using the original entries of the spectral density matrices and the sampling time At, as follows (Stengel, 1986): 941 = qlbt
+ qzAt3
(18)
1222 Ind. Eng. Chem. Res., Vol. 34, No. 4, 1995
(19)
rd = riAt
(21) The state space representation of the discrete filter will be
&&+I)= &(k)
- &)At
+ F(k)At + kdl[M,(k+l) A’&) + &@)At - F(k)At] (22)
+
k ( k + l ) = k(k) kd,[M,(k+l) - &(k)
+ &)At
-
F(k)Atl (23) The discrete transfer function in the complex variable z , which may be more appropriate for implementation, can be written as kd,[M,(Z)(Z
R ( z )=
z2 - (2 - kdl
-
1 ) - F(z)Atl
+ kdzAt)Z + ( 1 - kdl)
(24)
where kdl and kdz are related to the elements of the discrete state covariance matrix Pd = ( P d , } as in the COntinUOUS Case, i.t?., kdl =pdllhd and kdz =pdldrd. The discrete state covariance matrix at steady state can be obtained by solving the steady state discrete covariance equation (Stengel, 1986)
Pd
-
[I - pdHTRd-lH][@Pd@T+ Qd] = 0 (25)
where Q, is the state transition matrix of the discrete model, H is the measurement matrix, and &I and are the covariance matrices of the input and measurement noises whose entries are given by eqs 18-21. Although no analytic solution can be found for the covariance equation in this particular problem, after some manipulation it can be shown that the gains of the estimator defined above can be calculated after finding the roots of the following fourth order polynomial in a:
yy
At - a + A t 4(“:)2 - = O (26) The filter gains are related to the roots of this polynomial by
a2
k =I-dl
At2q4r kdz = a
(27) (28)
Note that the gains are already expressed as functions of the original spectral densities selected for the continuous estimator. Using Jury’s stability criterion (Astrom and Wittenmark, 1990) for the case when At(q2/r)1/4< 2, the only roots that give a stable filter meet the condition -At(q2/r)ll2 < a < 0. If the first inequality does not hold, another set of constraints on a can be obtained using the same criterion. 5. Applications
The estimation scheme will be applied to three examples of semicontinuous emulsion copolymerization.
In all cases the reactants (monomers) are added to the reactor a t a constant flow rate for about three-fourths of the total batch time. Then the feeds are stopped and the polymerization is allowed to proceed until most of the monomer is depleted. This type of operation is normally aimed to obtain a polymer of homogeneous composition. In order t o obtain a polymer of homogeneous composition,the rate of reaction of each monomer must be as constant as possible during the polymerization. This may be achieved by feeding each monomer at a very slow flow rate. In this case the rate of reaction of each monomer will coincide most of the time with the rate at which they are fed. If the feed rates are kept, low the chances of obtaining an homogeneous polymer increase; however low feed rates result in low productivity. In order to avoid such productivity losses, one wishes to operate at higher flow rates. For this reason it is desirable t o have a way of following the evolution of the reaction rates in order to make the system work as close as possible to the operating point that gives the highest productivity yet produce a homogeneous polymer. In the examples that will be presented no attempt to correct the evolution of the reaction rates is made and the reactions are run open loop. However the online monitoring of these rates is the first step toward the closed loop control of their evolution. 5.1. A Simulated Example. This example is based on simulated data from a very detailed model of the isothermal copolymerization of vinyl acetate (VAc)/butyl acrylate (BuA) (Dimitratos et al., 1991). This example will allow us to evaluate the estimator performance because the estimations can be compared with the “real” values. The model of the polymerization is used to generate 50 data points correspondingto measurements every 10 min of the total mass of VAc and BuA in the reactor. The flow rates of VAc and BuA are set at the beginning of the reaction to 0.946 x and 0.154 x lo-’ g/s respectively. After 380 min the flow rates are stopped and the polymerization allowed to proceed for an additional 110 min. The “measurements” are artificially corrupted with zero mean uniformly distributed random noise whose amplitude was chosen to be -5% of the maximum value of the noise-free “measurements”. Before implementing the estimator, the spectral density ( r )of the measurement noise and the acceptable variance ( p ) of the noise present in the estimation of the reaction rate must be established. The value of r is calculated according t o the selection of the artificial noise as explained above. For VAc r is 1.345 x g2 s whereas for BuA r is 2.3351 x lo-” g2 s. It is important to note here that the variance of the measurement noise (rd)and not its spectral density is what normally is obtained from calculations utilizing sampled data. To obtain the spectral density ( r )of the measurement noise, eq 21 must be used. In order to selectp for each reaction rate, the likely values of these rates must first be obtained. In this particular application (semicontinuous operation at low flow rates), the reaction rates are expected to be close in value to the flow rates of the corresponding components (monomers). Thus if a noise 5 10% of the likely value of the reaction rate is considered to be acceptable, the values ofp for VAc and BuA will be 2.237 x g2 s - ~and 5.929 x g2 s - ~ respectively. Using eq 14, values of qdr for VAc and BuA can be readily obtained. These values are qdr = 3.659 x sW4for VAc and q2/r = 2.986 x s - ~ for BuA. At this point it only remains to find the roots of the polynomial of eq 26. The polynomial must be
Ind. Eng. Chem. Res., Vol. 34, No. 4, 1995 1223 solved for each component. Out of the four roots for each one of the two components, one root will make the corresponding filter stable. Once these two values are obtained, the gains of the discrete estimators for both components can be readily calculated using eqs 27 and 28. For the first filter, the one that estimates the reaction rate of VAc, these gains are k d l = 0.4862 and k d z = -2.602 x sec-l. For the filter that estimates the reaction rate of BuA the pair of gains is k d l = 0.9774 and k d z = -1.557 x sec-’. Initial conditions are needed only for the state variables because the filter gains are fixed. Two initial conditions for each component i are needed. One i s for the estimated mass of monomer i in the reactor, Mj(O), and the other is forthe estimated rate of polymerization of that monomer, Ri(0). The selection of these values depends on how the polymerization is initiated. If no monomer is present a t time zero, both values should be set t o zero. In the case that a seeding stage is carried out before the semicontinuous operation, the initial conditions should be set according to this situation. In all our examples that comprised a seeding stage, the amounts of monomer were very small at the beginning of the semicontinuous operation. Consequently the polymerization rates were low. Therefore, the filter was always initiated with zeros in all states. The results of processing the simulated data with the filter described in eq 24 for the two pairs of gains according with the component being analyzed are shown in Figure 3. On the upper left hand corner of this figure the corrupted “measurement” of the mass of VAc in the reactor (+) is plotted together with the estimated value (Figure 3a). On the upper for the same variable (0) right hand corner of the same figure the “real” (-) and estimated (0) values of the reaction rate of VAc are plotted together with the flow rate profile used to feed VAc into the reactor (- - -1 (Figure 3b). Note that the noise in the estimated reaction rate is as predicted. It was explained in the theoretical section that the speed of the filter could be improved by allowing more noise in the estimation. Although the estimated value lags the “real” one at the ends, the overall result is quite satisfactory. For BuA the results are shown on the lower left and right hand corners, c and d, of Figure 3. The different variables are as described for the VAc plots. It is interesting to note that for the same amount of noise in the estimation the speed of the estimator for BuA is much faster than for VAc, producing an almost perfect estimate. In order to explain this behavior, let us assume for simplicity that the reaction rate to be estimated ( R )as well as the flow rate (F> are constant during the reaction time tf. Under this assumption
V(S)
= ~ W % (=Sy(F ) - R ) t ~ ( s ) (30)
where M”” is the maximum value of M during the reaction, y is the percentage of that maximum value that will determine the amplitude of the noise, and m(s) is the Laplace transform of a zero mean normalized random noise (e.g., a uniformly distributed random noise of amplitude 1). Let us now define
e =@2’/~
(31)
as the relative random error of the estimated reaction
x10.4
0
x10”
5
IO
TIME [HRS]
XlOd
C
0
5
TIME [HRSI
10
-5
-
0
5
10
TIME [HRSI
Figure 3. Theoreticalexample.(a) Mass of VAc (+) and estimated (b) Reaction rate of VAc (-1, estimated reaction mass of VAc (0).
rate of VAc (O), and flow rate of VAc (- - -). (c) Same as in (a) for BuA. (d) Same as in (b) for BuA.
rate, a value that should be kept small. Using the block diagram of Figure 2, it is easy to show that
e = ((FIR) - l>tfgL-’{G2(s) m(s)}
(32)
where L-l represents inverse Laplace transform and Gz(s) is the transfer function of the lower block in Figure 2. If for the same estimator we analyze two cases, one in which R = R1 and F = F1 and another in which R = R2 and F = F2, the ratio of the relative random errors will be given by
From Figure 2 it is clear that R1 and R2 will be tracked at the same speed because they are filtered by the same transfer function. However the relative amount of random noise in each estimation will depend on the ratio of eq 33. If F1IR1 > FdRz then el > e2, otherwise el < e2. This last finding helps explain the different behavior of the estimations for VAc and BuA. Using the argument derived assuming constant rates, the VAC/BuA example can be considered t o fulfill the relationship FVA~RVACFBUAIRBUA in which case eVAC > eBUA. Thus, the relative random error in the estimated value
1224 Ind. Eng. Chem. Res., Vol. 34, No. 4,1995
Figure 4. Schematic representation of the experimental emulsion polymerization reactor system: 1,jacketed glass reactor; 2,nitrogen supply; 3,injection port, 4,reflux condensor; 5, cold trap; 6, double impeller agitator; 7, baf’fles; 8, emulsifying charge inlet; 9, catalyst addition funnel; 10, catalyst addition pump; 11,jacket water heater and circulator; 12,internal cooling coil; 13,cooling water pump; 14, water buffer and overflow; 15,jacketed glass emulsifying tank; 16,impeller agitator; 17,cold water supply; 18,emulsifier addition pump.
would be higher for the VAc estimation if the estimators were designed to track the real values at the same speed. However, the amount of noise in the VAc estimation was selected at the design stage to be of the same order as that of the BuA estimation. Thus, keeping both relative noises at a similar level is achieved by sacrificing speed in the estimation of the reaction rate of VAc, as corroborated by Figure 3b,d. 5.2. Experimental Applications. The experimental applications were performed in a 1000 cm3jacketed glass laboratory reactor depicted in Figure 4 and also described in detail in (Dimitratos et al., 1991). The reactor has a bottom outlet for sample withdrawal and a glass four neck top section. To avoid uncontrolled polymerization, a nitrogen line was installed for blanketing the reaction with an inert atmosphere and removing the oxygen from the system during the reactor start-up. This line was also equipped with an injection port for initiator injection. A reflux condenser, with internal coil type cold finger was inserted into one of the side necks of the flask. The top of the condenser was connected to a one-piece cold trap, half-filled with water and open to the atmosphere for nitrogen pressure regulation in the reactor. The system was under atmospheric pressure. A double impeller agitator driven by a variable-speed motor (50-750 rpm) was employed to provide good mixing throughout the reactor. An open bath immersion circulator was used to raise the jacket
temperature and consequently the reacting mixture temperature to the desired value. The jacket acts as an insulator after the desired reactor temperature is reached. In order to control the reactor temperature, a stainless steel cooling coil was placed into the reactor. The flow of cold water through the coil was handled by a positive displacement gear pump. Six temperature measurements were taken at different locations in the reactor, the reactor jacket, and the coolant streams by using T-type thermocouples. Some of these measurements were used for closed loop control and some just for monitoring. Two inlet streams were used to provide the delayed ingredients during the reaction. One of the inlets was connected to a 50 cm3 addition funnel for initiator addition. This addition funnel has a nitrogen line installed to remove oxygen from the initiator solution. The other stream was connected to a emulsifying tank. This tank, which provides the emulsified monomers to the reacting mixture, is also a 1000 cm3jacketed glass flask with a glass four-neck top section. Cold water was circulated through the jacket of this flask in order to keep the monomer emulsion cool. A variable-speed stirrer provided the agitation needed to homogenize the naturally unstable monomer emulsions. A nitrogen line was connected to remove the oxygen from the monomer emulsions when needed. Two positive displacement metering pumps delivered the desired flow to the inlet
Ind. Eng. Chem. Res., Vol. 34, No. 4, 1995 1225 streams. The pumps were driven by computer-controlled stepping motors. Reactor temperature, jacket temperature, and flow rates of the additions and the cooling water were all controlled by an IBM compatible microcomputer using the process monitoring software ONSPEC. 5.2.1. Semibatch Copolymerization of Ethyl Acrylate and Acrylonitrile. In this application, the semibatch emulsion copolymerization of acrylonitrile (Acry)(15.19 partdl00 w/w) and ethyl acrylate (EA)(28 parts/lOO w/w) was conducted under starved conditions, that is with small amounts of monomer present during the polymerization, and isothermally at 83 "C. In addition to these two monomers, N-methylacrylamide in very small quantities (2.16 parts/100 w/w) was also used in the polymerization. Initially, part of the monomers (3.03 parts of Acry and 3.03 parts of EA), part of the initiator (0.33 parts of a 3% water solution of ammonium persulfate), part of one type of emulsifier (0.31 parts of a 30% water solution of sodium lauryl sulfate), part of another type of emulsifier (0.58 parts of a 35% water solution of Aerosol 22), and deionized water (29.34 parts) were polymerized in batch to produce a seed of polymer particles. After this seeding stage, the semicontinuous stage followed immediately. During the semicontinuous operation, the reactor was fed with (i) an emulsion of the remaining monomers in 0.55 parts of Sodium Lauryl Sulfate (30% water solution), 1.73 parts ofAerosol22 (35%water solution), and 16.72 parts of deionized water; and (ii) 5.42 parts of a 4% water solution of ammonium persulfate. After 3 h the flow rates of the monomer emulsion [Acry (0.00787 g/s) and EA (0.0163 g/s)l and the flow rate of initiator were stopped and the polymerization was allowed to proceed for another hour. During those 4 h, samples of the reacting mix were withdrawn every 10 min and analyzed using gas chromatography. A Hewlett Packard 5890 gas chromatograph equipped with a flame ionization detector (FID) and a fused silica column (0.53 mm i.d. x 10 m x 2.65 pm film thickness, cross-linked 5% PhMeSilicone, HP 5) was used for compositional analysis. A mass balance around the reactor and the use of butanol as internal standard allow us to translate the results from the chromatographic analysis into a quite precise quantitative measure of the total mass of Acry and EA in the reactor at every sampling time during the polymerization. The estimator in this example is designed to estimate the reaction rates of Acry and EA following the same steps described in the theoretical example. In this case the spectral densities of the measurement noises are calculated using data collected in previous experiments. First; the variances of a set of measurements of the mass of each component taken from a single sample are calculated for both monomers. Using eq 21, these values are transformed into spectral densities in order to finally design the appropriate estimators. The spectral densities ( r )of the noises are then calculated as 7.8 g2 s for the Acry measurement and 30.6 g2 s for the EA measurement. Using the same criterion as in the theoretical example, the variance p of the acceptable noise in the reaction rate estimates are obtained using the flow rates of each monomer. Therefore if a noise 5 10% of the likely value of the reaction rate is acceptg2 s - ~for Acry and 6.642 x able, p will be 1.548 x g2 for EA. With these values and using eq 14, qdr = 2.15 x s e c 4for Acry and qdr = 2.43 x ~ e c for - ~EA. Now solving the polynomial of eq 26 for
~10.3
TIME [HRSI
Y-
0
1
2
TIME (HRSI
TIME [HRSI
3
4
-
0
1
2
3
4
TIME [HRS]
Figure 5. AcryEA application. (a)Mass of Acry (+) and estimated (b) Estimated reaction rate of Acry (0) and flow mass of Acry (0). rate of Acry (-1. (c) Same as in (a) for EA. (d) Same as in (b) for EA.
a with the calculated values of qdr, the gains of the discrete filter for the estimation of the reaction rates sec-l for Acry are k d l = 0.970 and k d z = -1.523 x and k d l = 0.973 and k d z = 1.537 x sec-l for EA. The results of this experiment are shown in Figure 5. As in the previous example the upper two figures correspond to one of the monomers (Acry)and the lower two correspond to the other monomer (EA). The upper left part, a, of Figure 5 shows the measured amount of Acry in the reactor (+) and its estimated value (0). The upper right part, b, of Figure 5 shows the estimated reaction rate of Acry (0) and the flow rate of Acry that enters the reactor (-1. The same description applies for the lower figures except that they correspond to EA (Figures 4c,d). As can be observed, the reaction rates reach an almost constant value that coincides with the corresponding flow rates. These values are reached very quickly and are maintained until the feeding stage is over. After that, they decay rapidly following the flow rate profiles. Again the noise level in the estimations is as predicted by the theory. Although the real values are not available for comparison, these results are similar to the one obtained for the reaction rate of BuA in the simulated example. The small difference between the reaction rates and the flow rates anticipates an estimator with excellent tracking capabilities. 5.2.2. Semibatch Copolymerization of Butyl Acrylate and Styrene. This example is similar to the
1226 Ind. Eng. Chem. Res., Vol. 34, No. 4, 1995
previous one except that the monomers are butyl acrylate (BuA) (19.41 partdl00 w/w) and styrene (Sty) (17.47 partdl00 w/w). In addition to these two monomers, acrylic acid (AA)was used in small quantities (1.94 parts/100 w/w) in the polymerization. The temperature was kept at a constant value of 80 "C. Initially, part of the monomers (1.97 parts of BuA, 1.77 parts of EA, and 0.20 parts of AA),part of the initiator (0.01 parts of ammonium persulfate), part of the emulsifier (0.12 parts of a mixed surfactant), and deionized water (31.54 parts) were polymerized in batch to produce a seed of polymer particles. After this seeding stage the semicontinuous stage followed immediately. During the semicontinuous operation, the reactor was fed with (i) an emulsion of the remaining monomers in 1.08 parts of emulsifier and 28.39 parts of deionized water and (ii) 0.03 parts of ammonium persulfate; after 3 h the flow rates of the monomer emulsion [BuA (0.0133 g/s) and Sty (0.0120 g/s)l and the flow rate of initiator were stopped and the polymerization was allowed to proceed for another hour. The total amounts of BuA and Sty in the reactor were measured every 10 min as described in the ENAcry experiment. The variance of the measurement noises are calculated using samples from previous experiments. The corresponding spectral densities ( r ) are 2400 g2 s for BuA and 600 g2 s for Sty. The spectral densities ( p )of the noises in the estimations are selected as in the other two examples using the flow rates of both monomers as an estimate of the reaction rate values. If we assume that a noise of amplitude 5 10%of the value that is being estimated is acceptable, the variance p will be 4.42 x g2 s - ~for Sty. g2 s - ~for BuA and 3.6 x Following the same procedure as in the previous exs e c 4 for BuA and qdr = amples, qdr = 4.01 x 2.18 x sece4 for Sty and finally the gains of the discrete filter are k d l = 0.498 and k d 2 = 2.753 x lop4 sec-I for BuA and k d l = 0.644 and k d z -5.094 x sec-l for Sty. The results of this experiment are summarized in Figure 6. Measured and estimated values of the amounts of each monomer are shown in the left two parts, a and c, of Figure 6 and the estimated values of the reaction rates are shown in the right two figures together with the flow rate profiles of both monomers (Figure 6b,d). In this example contrary to what happened in the ENAcry application, the reaction rates approach very slowly the flow rate profiles. This result is consistent with lower values of monomer conversion found by gravimetric determination at the end of the reaction. The large difference between reaction rates and flow rate profiles anticipates a slower filter response as compared with the results of the ENAcry experiment. This result can be compared to that of the simulated example of the VAc monomer (see Figure 3b). In that example the feed and reaction rates of VAc are considerably different from each other during the polymerization. For reasons explained though the analysis that resulted into eq 33, the estimated rate lags in time behind the real one during the initial interval of the polymerization. During the remaining and major period of the polymerization, the filter performs its task quite accurately, as can be seen in Figure 3b. A similar performance of the estimator is expected to be taking place in the cases described by Figure 5b,d. Consequently one can argue that the proposed filtering scheme provides a satisfactory estimation of the reaction rates when the reactor is operated under conditions in
45
0.014
40
0.012 -
-
35
b
0.01 -
30
-
25
2 0.008-
20
0.006-
15
2 0.004-
E
x
10
TIME [HRSl
TIME [HRS]
1
0
1
2
TIME [HRS]
3
4
0
1
2
I
I
3
4
TIME [HRSI
Figure 6. BuA/Sty application. (a) Mass of BuA (+) and estimated mass of BuA (0).(b) Estimated reaction rate of BuA (0)and flow rate of BuA (-1~ ( c ) Same as in (a) for Sty. (d) Same as in (b) for sty.
which it is not possible to infer them directly from other operating variables, such as the correspondingfeed flow rates. 6. Conclusions
In this paper we have shown how linear Kalman filtering theory can be rigorously used to generate a useful tool to solve a practical problem arising in the monitoring of chemical reactors. The two degrees of freedom resulting from the application of the standard filter equations are used to our advantage to specify the desired filter behavior in simple terms. Out of the three parameters ( r , 91, and q 2 ) that typically appear in a design with one output and two states, the ratios ql/r and qdr were shown to govern the behavior of the estimator for a given chemical component. Because the input noise, w1,is fictitious and is introduced to model unknown variables, these two ratios are in principle adjustable parameters (Dalla Molle and Himmelblau, 1987). We have shown that using elements of control theory and stochastic systems theory they can be simply and uniquely determined. A major component of the present paper's contribution is the translation of the stability specifications found for the continuous estimator to a completely equivalent discrete filter design specification, which is accomplished through the solution of a fourth order polynomial.
Ind. Eng. Chem. Res., Vol. 34,No. 4,1995 1227 The theoretical results were applied first to a simulated example in order to test the ability of the estimator under different conditions. Then, the design methodology was applied to two experimental semibatch polymerizations. In both cases, the objective was to track the rates of reaction of two monomers which participated in an emulsion copolymerization of industrial interest. The ratio between these two rates determines the instantaneous composition of the produced copolymer, an important variable directly linked t o the properties of the final product. In all cases the results were satisfactory and the information needed to determine uniquely the estimators was easy to obtain.
Acknowledgment The authors gratefully acknowledge the financial support of the National Science Foundation under the IndustryNniversity Cooperative Research Centers Program and the financial support of Air Products and Chemicals as well as the other members of Industrial Consortium that supports research a t the Chemical Process Modeling and Control Research Center.
Dalla Molle, D. T.; Himmelblau, D. T. Fault Detection in a Single Stage Evaporator via Parameter Estimation Using a Kalman Filter. Znd. Eng. Chem. Res. 1987,26,2482-2489. Dimitratos, J.; Georgakis, C.; El-Aasser, M. S.; Klein, A. An Experimental Study of Adaptive Kalman Filtering in Emulsion Copolymerization. Chem. Eng. Sci. 1991,46(121,3203. Lewis, F. I. Optimal Estimation; John Wiley: New York, 1989. Min, K. W.; Ray, W. H. On the Mathematical Modeling of Emulsion Polymerization Reactors. J . Macromol. Sci., Rev. Macromol. Chem. 1974,C11 (2),177. Papoulis, A. Probability, Random Variables, and Stochastic Processes; McGraw-Hill: New York, 1984. Ray, W. H.; Gall, C. E. The Control of Copolymer Composition Distributions in Batch and Tubular Reactors. Macromolecules
1969,2(41,425. Schuler, H.;Schmidt, C.-U. Calorimetric-State Estimators for Chemical Reactor Diagnosis and Control: Review of Methods and Applications. Chem. Eng. Sci. 1992,47, 899-915. Stengel, R. F. Stochastic Optimal Control; John Wiley & Sons: New York, 1986. Tirrell, M.; Gromley, K. Composition Control of Batch Copolymerization Reactors. Chem. Eng. Sci. 1981,36,367.
Received for review May 31,1994 Revised manuscript received December 2, 1994 Accepted January 10,1995@
I39403479
Literature Cited Astrom, K.; Wittenmark, B. Computer-Controlled Systems. Theory and Design; Prentice Hall: Engelwood Cliffs, NJ, 1990.
@
Abstract published in Advance A C S Abstracts, March 1,
1995.
