I n d . Eng. Chem. Res. 1989,28, 1846-1861
1846
Continuous Solution Polymerization Reactor Control. 2. Estimation and Nonlinear Reference Control during Methyl Methacrylate Polymerization Derinola K. Adebekun and F. Joseph Schork* School of Chemical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0100
This paper focuses on the estimation and control of a continuous reactor during solution polymerization of methyl methacrylate in a CSTR. The leading moments of the molecular weight distribution are estimated. Various alternatives for the control are considered, and the effects of poor filter initialization are simulated. T h e results indicate that these estimates could be adequate for feedback control of the polymerization reactor.
Introduction Polymeric materials are employed in a vast number of applications. The end use of these materials is to a large extent dependent on the molecular weight distribution (MWD) of the polymer. This distribution can be controlled only during the manufacture of the polymer and is a complicated function of the reactor conditions. Hence, good control of the MWD motivates the need to monitor the reactor state accurately. Usually, for quality control purposes, the MWD is characterized by the number-average chain length ( w ) and the polydispersity (D). These variables can be inferred from the leading moments of the distribution (Ray, 1972). Typically, not all the state variables of the reactor required for control are measurable. Even when they are, such measurements might be infrequent. Thus, there is a need to reconstruct unmeasurable (and/or infrequently measurable) state variables from the available measurements. To this end, techniques in optimal estimation are finding increased applications in polymerization reaction engineering in the form of the extended Kalman filter (EKF). The first application of Kalman filtering in polymerization appears to be the work of Jo and Bankoff (1976). In that study, good estimates of conversion and average molecular weight were obtained from refractive index and viscosity measurements during CSTR polymerization. More recently, attempts have been made to estimate the entire MWD during batch polymerizations. Recent advances in the area include the work of Ellis et d. (1988a,b), Schuler and Suzhen (1985),and Papadopoulou and Gilles (1986). In particular, in some of these works, the authors incorporate time-lagged gel permeation chromatography (GPC) measurements. This has led naturally to the development of the two-time-scale filtering algorithm (Ellis et al., 1988a,b) and the application of some model reduction techniques. Other related work on Kalman filtering has been presented by various authors (Choi and Khan, 1988; Dimitratos et al., 1988, 1989). In this work, we apply some of these techniques to solution polymerization of methyl methacrylate (MMA) in a CSTR during open-loop dynamic simulations. We demonstrate the viability of these estimators for this system as well as the robustness of the filter to poor initialization. In so doing, we also show that, even in the event of incomplete state observability, the open-loop performance of these estimators could still be acceptable.
* Author
to whom all correspondence should be directed.
We go one step further by employing these estimates in closed-loop control schemes as very few such studies are available. A notable example involving linear optimal control theory is the paper by Kipparissides et al. (1981). In the first paper (Adebekun and Schork, 1989a) of this two-part series, henceforth referred to as part 1 , we developed and analyzed the convergence properties of certain nonlinear reference controllers designed for the deterministic plant. These controllers required the availability of the state variables of reactor subsystem 1 (consisting of monomer, temperature, initiator, and solvent levels) for implementation. Thus, we estimate the required state variables and use these for control even in cases where the desired equilibrium point is open-loop unstable.
Model Methyl methacrylate (MMA) is commonly polymerized via a free-radical mechanism. In this study, the solvent is ethyl acetate while the initiator is benzoyl peroxide. The model assumed is, in essence, taken to be that of Schmidt and Ray (1981). As in Tanner et al. (1987), constant density is assumed and the equation for the second moment of the MWD is included. As was discussed in part 1, the nonlinearities in this model are quite severe. We note that the model includes terms for chain transfer to both monomer and solvent. As is well-known, chain-terminating reactions should be considered when the MWD is of interest. Also, we further observe that for MMA, at least, arguments can be made for including terms that take into account pure thermal initiation (Odian, 1981). We, however, do not include this term in the spirit that, in practice, the EKF might capture the effects of such terms via process noise inclusion. Indeed, as Schuler and Suzhen (1985) demonstrate, certain terms could be omitted in the model with negligible degradation of filter performance. Additionally, even though termination is mainly by disproportionation, transfer to polymer can be safely ignored for MMA (Odian, 1981). Thus, as in part 1,the following equations can readily be derived:
VM = q ( M f - M ) - VkJklP
(1)
VpcPT = pcpq(Tf- T ) + (-AH)VkJklP - hA,(T - T,) (2)
V i = q(Zf - I ) - VkdZ
(3)
vs = q ( s f - S )
(4)
0888-5885/89/2628- 1846$01.50/0 0 1989 American Chemical Society
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1847 Vi0
= -qX,
1 + V[(kfM + k,dP + kf8S)]~P + iVk,cP2
(5)
vi1 = -qX1 + V[(kfM + kaP + kf8S)(2a - a2)+ k,cP]P/(l - a) (6) Vi2 = -qhz + V[(kfM+ kaP + kf,S)(a3 - 3a2 + 4a) + k,cP(a + 2)]P/(1 - a)2 (7)
Dap = k,'elPVMfo/q
Daa = ka:e-Y@pVMfQ/q
kdM
+ kfM + kf,S + ktP
(0.10575 exp[17.15Vf - 0.01715(T - 273.2)]
V, > [0.1856 - 2.965
X
10-'(T
- 273.2)]
(9)
2.3 X lo4 exp(75Vr)
\ Vr I [0.1856 - 2.965
X
10-"(T - 273.211
where Vf is the free volume of the polymer calculated from the volume fractions of monomer, polymer, and solvent in the reactor (Schmidt and Ray, 1981). Thus, as in Jaisinghani and Ray (1977),the following dimensionless variables are introduced:
+ BDa,y,Wx,E,
- p(xz - xZc)
13 = X3f - X3 - DUdX3Exd
(8)
P = (2flz,J/k,)1/2 is the total concentration of live radicals, and k,(=k, + ktd) is the overall termination rate constant. The ith moment of the dead polymer MWD is represented by X i (i = 0, 1,2), and M is the monomer concentration. T, I, and S are the reactor temperature, initiator concentration, and solvent concentration, respectively. The rate constants with the exception of kf, are assumed to follow an Arrhenius dependence on temperature. We emphasize that, quite apart from the nonisothermal effect, the gel effect introduces severe nonlinearities into the model equations. As discussed in part 1, it is directly responsible for the occurrence of input multiplicity during closed-loop control for this reactor, and this can lead to difficult control problems in practice. Furthermore, it leads to complicated open-loop steady-state (output) multiplicity patterns (in the classical sense, as distinct from input multiplicity) in these reactors (Schmidt et al., 1984; Adebekun et al., 1989). Thus, any realistic attempt to model this reactor should include at the very least, a term to account for its presence. To this end, we employ the Schmidt-Ray (1981) correlation for the gel effect. The relevant equations are
Da, = k;e-YflPVMfo/q
and the dimensionless model equations become xl = xlf - x1 - Da,Wx,E, 3t2 = -x2
k$4
Da, = kfBMfQV/q
Da, = ktc:eYkypVMfQ/q
Da, = kt:e-YtYPVMfQ/q
where a =
Dad = vDaP
14
3a2 + 4a)
= X4f - x4
+ g,Da,WE,(a
(10) (11)
(12)
(13)
+ 2)] W/(l
- a ) 2 (16)
where
and
Some comments are in order concerning the glass effect correlation (g,) sometimes employed for the propagation rate constant. We set this factor to unity in this work as this is often known to be a good assumption for moderate solvent concentrations. The model above represented by eq 10-18 is thus believed to be a good representation for the system and represents the nominal plant chosen for this investigation. The model parameters and kinetic data are shown in Table I.
Results and Discussion The results are organized as follows: in section I, the control schemes are briefly described and further clarifications are made with regard to the results obtained in part 1. Section I1 contains the results on some system theoretic properties of the system. This is done with the hope of providing some theoretical basis for the simulated estimation algorithms. In section 111.1, we present the first estimation approach and show some simulation results. The algorithm in this section is referred to as a full-order estimator in a sense to be made precise later. In section 111.2, we modify the algorithm in section 111.1 to obtain, in essence, the Schuler-Suzhen (1985) approach in which we identify subsystems 1 and 2 of the reactor explicitly. Section 111.3 addresses the information injection problem of introducing off-line MWD measurements into the process. Note that, throughout section 111, we focus on the open-loop plant only. Finally, in section IV, we direct our attention toward interconnecting the state feedback developed in section I and the state estimator. I. Deterministic Control. 1.1. Least-Squares Control. We consider the problem of driving the reactor from
1848 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
Table I. Datao lit. values (mol, L, s) ref Kinetic Rate Constants kd 1.69 X lo1*exp(-30000/RT) 4.925 X lo6 exp(-4353/RT) k, 9.80 X 10' exp(-701/RT) k, kf 4.92 exp(-4353/RT) Brandrup and Immergut, 1975 kfs 0.091 kdlk, 8.23 Reaction Medium Parameters Schmidt and Ray, 1981 f = 0.5 h = 135.6 cal/(m2.s.K) A, = 2.8 m2 p = 1038 g/L Rodriguez, 1982 Jaisinghani and Ray, 1977 c p = 0.4 cal/(g.K) Jaisinghani and Ray, 1977 -AH = 13.8 kcal/mol kj
"Table I: Schmidt and Ray (1981) unless otherwise noted.
some initial condition (x,) to some open-loop equilibrium point (xd). As usual, the manipulated variables (u E R 3, are chosen to be the monomer inlet concentration (xlf), the jacket temperature (xzc), and the initiator inlet concentration ( x d . Observe that the model equations are of the form x = f(x) Gu (19)
+
where
four components. As was discussed in part 1, we will switch to a seven-component state vector when we consider the MWD explicitly. The results on the stability and convergence of the system to xd as well as the choice of an appropriate k are available in part 1. 1.2. Exact Linearization of a Subset of the State. Here, we focus on controlling a subset of the state, specifically the first three state variables. Even so, as was proven in part 1,under certain conditions, global stabilization of the entire plant can also be achieved. As before, define a reference system for the first three state variables (monomer, temperature, and initiator levels) of subsystem 1: 211= Al(x11- xdl)
(24)
where A, E R 3 X 3and xI1and xdl E R 3 . The vector xdl is the portion of the desired equilibrium point consisting of the desired levels of monomer, temperature, and initiator, respectively. Again, matrix Al is chosen to be diagonal and of the form -kI with k > 0. Then, the dynamics of the error system e , (=xl - x,,, E R 3 ) are given by e, = Alel + [f,(x) - Al(x, - xdl) + Glu] (25) The vector x1 E R defines the first three state variables (monomer, temperature, and initiator levels) of the reactor. fl(x)is the vector consisting of the first three rows of f(x) in eq 19. (Remark: fl(x)depends on all the state variables of subsystem 1-dependency on x4 comes in through the gel effect factor gt). The nonsingular matrix G1 (E R 3x3) consists of the first three rows of the matrix Go in eq 19. Again for e 0 as t a,we equate the term in the square bracket in eq 25 to zero so that we have the exact solution:
- -
The vector f(x) E R 7 in eq 19 represents the remaining terms in eq 10-17 after the above rearrangement. Material balance considerations ensure that a, and a3 are nonzero since the solvent feed concentration (xlf) is not independent of xlf and xQf. Also, /3 # 0 and the matrix GTG is nonsingular. A first-order response is desired in all loops. Thus, we define a reference system as X r = A(X, - xd) (20) where the matrix A is chosen to be diagonal, asymptotically stable, and of the form -kI with k > 0. By defining the error system e = x - x,, we note from eq 19 and 20 that e = Ae [f(x) - A(x - xd) + Gu] (21)
+
- -
In ensuring that e 0 as t a,we need to determine u such that O = f(X) - A ( X - xd) + G U (22) The dimensions of the various vectors defined earlier (more equations than unknowns) show that we reduce the problem to a least-squares estimation problem so that u = -(GTG)-'GT[f(x)- A ( X - xd)]
(23)
Equation 23, which is defined for all times, is a leastsquares solution to eq 22. Remark. The structure of G shows that the state variables of subsystem 2 (the MWD subsystem) do not appear in the control law. In essence, we could have obtained the same control law by merely considering subsystem 1 alone and defining a corresponding four-state reference system. Thus, in the rest of this paper, we simply replace (in a Fortran sense) G by Go by assuming G E R 4x3. This assumes that the state vector (x) consists of
U
= -Gl-'[fi(x) - Al(x1 - XdJ]
(26)
In part 1, the above control law is proven to be stable (if k is chosen to ensure that the physically unrealistic condition Q I 0 is not encountered in closed loop; see section 1.2.1 below). However, using this scheme, the closed-loop Jacobian shows that it may be impossible to stabilize certain members of the catastrophe set. Furthermore, for all such k, global stabilization is possible for all equilibrium points not belonging to the catastrophe set. If the condition Q > 0 is violated, it is possible in theory for the inputs to converge to physically unrealistic values, a situation that in practice implies instability. Indeed, as is often done in the literature, if the gel effect (gJ is neglected, the catastrophe set is nonexistent and asymptotic stability is global for all equilibrium points (V k ) . (Remark: Q I0 is a major, but not the only, reason for convergence to unrealistic inputs.) 1.2.1. Practical Considerations in Control. At this point, it is necessary to further clarify some details concerning the exact linearization controller as they relate to practical considerations. These comments also address some of the input-multiplicity issues raised in part 1. The stability result in part 1 (section 11.1) depends crucially on the fact that the free volume (Vf) is bounded below by zero. In theory, we accomplished this by a straightforward extension of the gel effect correlation of Schmidt and Ray (1981) in which we introduced an auxilliary variable (Q)defined by (27) = Vfnl&ll + vfp4p + Vf.qd4 so that Vf =
Q , if Q > 0 0, if Q 5 0
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1849
.o
n l
0.0
3.0
i.0
i.0
*.O
li.0
14.0
I
8.0
9.0
(2.0
15.0
8.0
8.0
11.0
(5.0
n
E % n
0.0
3.0
Resldence time unlts
Figure 1. Full-order estimator: entire state vector estimated from measurements of conversion and temperature; 1denotes true plant while 2 represents the estimate.
Hence, V, exactly captures physical reality in the region R > 0. The second condition in eq 28 in effect implies that we assumed a minimum nonzero value (gt,min)for the gel effect correlation, reasonable since the chains have to terminate. Furthermore, in order to exclude the input corresponding to g, 1, we also assumed g , < 1 for R > 0. In effect, the correlation is always active after start-up and at equilibrium. (Remark: Define Vf 0.025 if R 5 0 does not qualitatively affect the following discussion.) However, the second condition (in eq 28) always introduces in theory, a unique spurious extraneous steady-state value (x4,min)for the x4 (solvent) component for several equilibrium points when the first three components are maintained at the desired values, xdl. This is because, in the region Q I 0, we have identically that g , gt,min.Because gmh is (at most) on the order of 2.3 X lo+, it is not too difficult to argue that this spurious equilibrium point is almost always not in physical parameter space (certainly not at the conditions simulated in case study 6-Figure 13 of part 1). This is why we were able to conclude in that case study that the system was seeking for this equilibrium point. Furthermore, because this equilibrium point is always stable, if D < 0, it could be encountered in closed-loop dynamic simulations (Figure 13, part 1)if the controller gain (k)is such that, at some point, it drives the closed-loop plant into the region Q I 0. In static bifurcations, when we deal with the open-loop plant, analysis is restricted to the region of model validity ( Q > 01, in which case this spurious equilibrium point is never encountered and, thus, should not be included in considerations involving the catastrophe set. In effect, for practical purposes, the catastrophe set is defined for inputs only when > 0; a degeneracy occurs
(as in Figure 13, part 1) if the region D I 0 is allowed. In addition (ignoring other physical constraints), k has to be chosen to ensure that this restriction is always satisfied. All these points, even though implicit in part 1, should have been explicitly stated in order to enhance clarity. The above discussion then clearly indicates the limitations on the range of acceptable gains @)-this was why we usually picked k to be low in part 1. In practice, one would choose k low enough to ensure that a small neighborhood of the condition Q = 0 is not encountered in closed loop. As the simulation results in part 1show, this is not a serious problem. For such k values, then, global stabilization is possible for equilibrium points not belonging to the catastrophe set. Indeed, physically, a violation of this condition implies that the model has broken down; in effect, in practice, the closed-loop system has become unstable. It is interesting to note that this problem arises because we exit physical parameter space-due to the gel effect correlation. It does not arise in several other reaction schemes such as in nth-order nonisothermal reactions ( n E (0,m)). Indeed, it is proven that any equilibrium point of such reactions belonging to physical parameter space is globally stabilized by this class of controllers (Adebekun and Schork, 1989b). 11. System Theoretic Properties. In any estimation scheme, it is important to consider the observability properties of the system. This need arises from the fact that most of the theorems concerning the convergence properties of the filter demand system observability (for linear systems). Actually, it is often possible to devise good estimation algorithms if the less restrictive detectability requirements of the system are met. The difficulty is that observability is more easily verified than detectability. A
1850 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
-E
E
9-
n-
.O
Rssldence t l m s u n i t s
Figure 2. Full-order filter: effect of process noise increase; 1, true plant; 2, estimate.
1/ '
.o
1.0
0 2.0
3.0
4.0
5.0
8.0
I
7.0
0
0
x ;
c 1 e
-n
*-
/ a ? a L 0-
p.---
EE
--
I L
L
E:-
$2
-E d 2" c . n I
'b: f i-A/ e
L ?
E n
eo
E?
E m m
- 0
a a
,
,
/Hl
5; E
:-
:s -,
w 0 c -0 o-
m0.
-
C O
zi-l/
0 -
a
0.0
1.0
1.0
3.0
1.0
5.0
1.0
7.0
-aE 0d ,
Residence t l m s unlts
Figure 3. Reduced-order filter with perfect knowledge of P; 1, plant; 2 estimate.
/LL ,/
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1851
I
.
.
!
.
.o
J.0
i.0
i.0
12.0
0.0
i.0
i.0
i. o
12.0
I
15.0
t -
:i
-* C0 . -
n
5:
I,
--.---.,
,
,
,
,
,
,
,
,
. 0
0.0
1.0
, 1.0
8.0
,
, lZ.0
,
11.0
Rerldrncr tlmr unltr
Figure 4. Reduced-order filter, -5% change in 8; 1, plant; 2, estimate.
very good discussion and examples of applicability of these theorems are presented by Lewis (1986). Indeed, for multivariable linear time-invariant systems of the form x = Fx Bu G'w (29)
+
+
y=Hx+v (30) where w and v are normally distributed, zero mean, white noise processes, the convergence of the corresponding continuous-continuous Kalman filter is guaranteed provided R > 0, (F,H) is detectable and (F,G'Q1/2)is reachable. In eq 29 and 30, vectors x, u, and y represent the system state, deterministic input, and system outputs, respectively. Matrices R and Q denote the covariance of v and spectral density of w, respectively (Jazwinski, 1970; Lewis, 1986;Anderson and Moore, 1979). In order to avoid trivialties, all matrices are assumed to be appropriately dimensioned. The above result provides some justification for some of the designs which follow. Rigorously, because the system is nonlinear (hence, the Jacobian is time varying) and an EKF is being employed, these results are not exactly applicable to our system. However, as will be later seen, the (Taylor) linearized system meets the above requirements. In particular, the structure of the Jacobian makes the verification of the detectability condition trivial. Observe from the model equations (10)-(17) that the linearized system has a Jacobian F of the form
where Ali E R4X4and Az2E R3x3. Matrix Azl E R 3 X 4 and the dimensions of the 0 block can easily be inferred. The important point to note is that F is block-triangular. Thus, we reduce to the fact that the eigenvalues of F are
equal to those of matrices All and AZ2(Owens, 1978). Furthermore, it is clear that all the three eigenvaluesof AZ2are equal to -1. Thus, as is quite well-known (Ray, 1986), the MWD modes are aymptotically stable. Next, we consider matrix All to see that this is of the form A1111 = (0
Alll* A,,,,)
(32)
which is also a block-triangular matrix. Hence, its eigenvalues are given by those of Allll (E R3x3)and AIlz2(E R). It is easy to see that All = -1, in which case the mode corresponding to the solventloop (x4) is also asymptotically stable. Our methods of state estimation will depend on the above reduction. Remarks are in order here about filter tuning. The initial tuning guidelines were always based on choosing the noise matrices to satisfy the detectability and reachability conditions discussed earlier. Additional fine-tuning was sometimesachieved by performing numerical experiments. These issues are further discussed as they arise. 111. Approaches to the State Estimation Problem. In this section, various approaches to the state estimation problem are implemented. Three schemes are considered, and these are described as they arise. In the first approach, we estimate the entire state vector (with seven components) by using available measurements. For obvious reasons, we refer to this implementation as a full-order estimator. The second method involves some model reduction similar to the method of Schuler and Suzhen (1985), in which we estimate the subsystem 1state vector and do a pure prediction on subsystem 2. In the third approach, we implement the two-time-scale filter which considers the situation when, in addition, delayed measurements of subsystem 2 are available.
1852
.o
0
E
f
5-
$-,
=-.
c . I?.
-E o
/i,,
E I
E?
- 0
q+ii, 4
-a:€::2
0.0
1.0
2.0
1.0
4.0
5.0
6.0
7.0
8.0
s.0
10.0
.
