Ind. Eng. Chem. Res. 1999, 38, 4815-4824
4815
Applications of Nonlinear Observers and Control: Improving Productivity and Control of Free Radical Solution Copolymerization H. Hammouri,*,† T. F. McKenna,‡ and S. Othman† LAGEP-UPRESA CNRS Q5007/CPE, University Lyon I, 43 Blvd. du 11 Novembre 1918, 69622 Villeurbanne Cedex, France, and LCPP-CNRS/CPE, Baˆ t. F308, 43 Blvd. du 11 Novembre 1918, 69616 Villeurbanne Cedex, France
A nonlinear state observer that uses the rate of heat generation due to chemical reaction for the real-time evaluation of key parameters during free radical solution copolymerization is presented. These parameters include the cumulative and instantaneous copolymer compositions, the overall and individual monomer conversions, and the overall termination rate constant associated with the reaction. The information thus obtained is used to control the cumulative copolymer composition and to maintain the productivity of a semibatch reactor for the solution copolymerization of butyl acrylate (BuA) and vinyl acetate (VAc) in ethyl acetate (EAc) at its highest possible level. Both the state and parameter estimation of this observer were then used for the nonlinear control of copolymer composition and optimal reactor productivity. 1. Introduction and Objectives Typically, the task of a process engineer is to simultaneously maintain product quality and high production rates. While this can be a difficult task for any process, it is particularly difficult in the case of polymerization reactions, given the highly nonlinear nature of the relationship between concentrations and reaction rates, the difficulty of monitoring properties related to enduse (e.g. composition and molecular weight) on-line, and the fact that it is very difficult to obtain precise values of complex rate constants that are needed to model the system in question. With the current trend toward more highly computerized systems, it is useful to take advantage of this cheap, rapid computation21 power now available in developing solutions to such problems. Such efforts involve the combination of hardware and software sensors that are connected to the polymerization reactor and, with the use of more or less precise models of the system in question, allow us to obtain information on the evolution of the reaction rate and so forth when combined with powerful, robust state estimators and control laws. We will not discuss this point in more detail here, as it has been treated in a large number of studies over the course of the past few years. For further information the reader is referred to refs 1-5 (while this list is far from exhaustive, it will provide the reader with a good overview of the current state of the art). The control of polymerization reactors has been the object of a great number of publications. A relatively recent article6 presents a review of polymer reactor control, as well as of certain on-line measurement techniques including calorimetry, and a different review article covers some of the key issues in polymer reactor control that include the need for realistic physical models.7 There are essentially two broad categories of approaches to the construction of nonlinear observers (or estimators). Most of the control methods discussed * To whom correspondence should be addressed. E-mail:
[email protected]. † LAGEP-UPRESA CNRS Q5007/CPE. ‡ LCPP-CNRS/CPE.
in this article involved extended Kalman filtering.8-10 The Kalman filter-based approach is based on a linear approximation of a nonlinear process model. This method can produce satisfactory process control schemes, but in most cases, there is no a priori guarantee of the convergence and stability of these algorithms.8,11 The estimators proposed in earlier works8-10 seem to present good convergence properties and robustness to both state and measurement noise. However, estimator tuning is not a trivial exercise. On the other hand, the nonlinear estimation and control techniques presented here have the advantage of being relatively simple to tune. As was discussed elsewhere,3,4 there are far fewer tuning parameters with this nonlinear approach than with the Kalman filter. Other nonlinear approaches have also been discussed in the literature.12 In this work, Kravaris et al. used nonlinear temperature tracking to control the copolymer composition batch and semibatch processes. Although their approach appears to be useful for maintaining constant composition in end-of-batch finishing periods, its open-loop state estimation requires that one be entirely certain of the proces modelssomething that is not absolutely requried of the approach here (however, the better the model, the more informaiton we can obtain from the process). In addition, they limit themselves to composition control. The work we discuss below shows how to use high-gain observers to simultaneously control the composition and production rate. In so far as calorimetry is concerned, it has been demonstrated by various research groups that this is a very powerful means of extracting information from a polymerization reaction.2-4,13-19 However, it has been pointed out3,4 that the problem with a large number of the approaches in the literature is that one needs to understand how the overall heat-transfer coefficient varies as a function of time, as well as to make assumptions about heat losses and so forth. This is not true of the calorimetric approach used in this work. This method, based on a cascade of two state estimatorss one that is used to close the energy balance and one that provides information about the state of the reactor
10.1021/ie9806996 CCC: $18.00 © 1999 American Chemical Society Published on Web 11/05/1999
4816
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
so we will simply present the model equations used here. We will consider only semibatch systems in the current work, as they are common when composition control in multicomponent systems is important. As in most engineering applications, a description of the system relies on mass and energy balances. For a semibatch reactor such as the one shown in Figure 1, the evolution of the number of moles of each species, that is, monomer (N1 and N2), initiator (NI), and solvent (NS), and the volume of the reactor contents (if we ignore the small contribution of the change of volume on reaction) are given by Figure 1. Reactor scheme showing two separate, independent control loops: one that maintains the isothermal reactor temperature and one that controls the flow rates of the monomer and the initiator.
dN1 ) q1,in - Rp1VP ) q1,in dt [R*]N1(kp11φ1* + kp21φ2*)
contentsshas been tested experimentally in solution4 and emulsion systems.20,21 In the current work, we will address the specific case of the free radical solution polymerization of two monomers. The objectives of the control scheme discussed here are two-fold: one, to maintain the cumulative composition of a copolymer at a specified value (which might or might not vary as a function of time) and, two, to improve reactor productivity. Such a control scheme will necessarily be based on a robust state observer that will allow us to know with certainty what is happening in the reactor at any given time. It was clearly shown4,20,21 that this method is valid in the presence of inhibitors, nonideal kinetic behavior, and nonisothermal reactor operation. For this reason, we will use simplified energy balance equations and kinetic expressions in the current work to demonstrate the applicability of the control laws, knowing that we could obtain the necessary information from a more complex system should the need arise. Temperature control of the reactor is important, since free radical polymerizations are highly exothermic. On one hand, high rates of heat release make calorimetry a useful monitoring tool. On the other hand, the rates of reaction, and thus copolymer composition and molecular weights, are highly temperature dependent. This, combined with obvious safety concerns, means that heat removal is an important task. The maximum amount of heat that can be removed from a chemical reactor is finite and can generally be estimated rather easily. If one wishes to maintain isothermal conditions (or to avoid thermal runaway of the reaction), this limiting heat removal rate must not be exceeded. However, for semibatch systems it is not uncommon for the heat production rate, and therefore reactor productivity, to be below their maximal values for long periods during the reaction.22 To overcome this weakness, one can continuously feed controlled quantities of initiator to the reaction as a function of time in order to maintain the heat production rate and thus productivity at the highest acceptable levels.22 In the current work, we will combine both of these aspects (a robust state observer and the maximum productivity based on heat removal capacity) in order to provide targets for our nonlinear control scheme.
dN2 ) q2,in - Rp2VP ) q2,in dt [R*]N2(kp12φ1* + kp22φ2*)
2. Kinetics and On-Line Measurement 2.1. Kinetic Model. The kinetic scheme for free radical polymerization and the basic mass and energy balance equations have been described elsewhere,23,24
dNI ) qI,in - kdNI dt dNS ) qS,in dt dVp
)
dt
∑i
qi,inMWi Fi
(1)
i ) 1, 2, I, S
where q1,in and q2,in are the inlet molar flow rates of monomer. Similarly, qI,in and qS,in are the inlet flow rates of initiator and of solvent, kp11 and kp22 are the homopolymerization propagation rate constants for monomer 1 and monomer 2, and kpij is the cross-propagation rate constant between a radical of type i and monomer j (i * j). The cross-propagation rate constants are generally estimated from the homopolymerization constants divided by the well-known reactivity ratios. In the case of a copolymerization, these ratios are defined as
r1 )
kp11 kp22 and r2 ) kp12 kp21
and φ1* and φ2* are the mole fractions of radicals terminating in by a radical of monomers 1 and 2:
(
φ1* ) 1 +
)
kp12 N2 kp21 N1
-1
and φ2* ) 1 - φ1*
(2)
[R*] is the total concentration of free radicals (we have of course implicitly neglected the contribution of initiator radicals). If the quasi-steady-state assumption is valid, and in the absence of inhibition, the total radical concentration for a single initiator is given by
[R*] )
( ) 2fkd[I] kt
1/2
(3)
Here, [I] ) NI/Vp is the initiator concentration, kd is the initiator decomposition constant, and, finally, kt is the overall termination rate constant:
kt ) kt,11φ1*2 + kt,22φ2*2 + 2ωxkt,11kt,22φ1*φ2* (4)
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4817
kt11 and kt22 are the homopolymerization rate constants that can be found in the literature. Note that the expression of kt is widely used in the literature, although there is very little evidence that ω can be taken as a constant. For this reason, coupled with the fact that the values of the homopolymerization rate constants vary as a function of polymer concentration and molecular weight in rather unpredictable ways, it is very useful to estimate kt as it evolves over the course of the reaction. 2.2. Calorimetry. The use of calorimetry as an online technique in polymer reaction engineering has been discussed in the literature for several years now.3,4 The operating principle is simple: an energy balance is performed around the reactor, and all known or measured terms are used to estimate the heat of reaction. The energy balance is usually composed of the following terms:
Qacc ) mrcpr
dT ) QF - Qcond + QR + Qm - Qj - Qloss dt (5)
where Qacc is the rate of accumulation of energy in the reactor, mrcpr is the product of the mass of the reactor contents and an average value of the heat capacity (this includes any reactor fixtures such as the agitator or baffles), QF is the sensible heat of the reactor feed, Qcond is the heat removed in the condenser, Qm is the mechanical energy input due to agitation, Qj is the heat removed through the cooling jacket, and Qloss is the term representing the heat lost to the surrounding environment. The heat of reaction is related to the rates of polymerization, and ultimately the individual conversions of monomer, by
QR ) Vp
∑i
2
∑ i)1
Rpi(-∆Hpi) ) [R*]
2
{-∆Hpi
kpjiφj*Ni} ∑ j)1
(6)
This expression shows the explicit relationship between the radical concentration, the monomer composition, and the heat of reaction. Since we have given ourselves the task of composition control, the first question is therefore how to define copolymer composition. There are in fact two different types of composition: the instantaneous and the cumulative compositions. The instantaneous composition refers to the composition of the copolymer produced at any given instant in time (the relative rates of entry of monomers 1 and 2 into the polymer) and is the ratio of the rates of reaction:
{
Rp1 r1f12 + f1(1 - f1) F1 ) ) Rp1 + Rp2 r f 2 + 2f (1 - f ) + r (1 - f )2 1 1 1 1 2 1 F2 ) 1 - F1 (7)
where f1 is the mole fraction of monomer 1 relative to monomer 1 plus monomer 2 f1 ) N1/(N1 + N2) and r1 and r2 are the reactivity ratios defined above. This relationship between F1 and f1 is important. In fact, in a typical example of composition control, one specifies F1(t) (usually constant) and needs to obtain a specific value of f1. The control policy is therefore to act on maintaining a set profile for f1 as a function of time (e.g. constant if F1 is constant). Finally, if we can estimate
the variation of the number of moles of all of the monomers in the system, the cumulative copolymer compositions are calculated from the following expressions:
F1 )
x1N1tot x1N1tot + x2N2tot
(8)
where xi is the individual monomer conversions of monomer i given by xi ) Nitot - Ni/Nitot, and Nitot is the total number of moles of monomer i added to the reactor Nitot ) Ni(0) + ∫t0 qi(τ) dτ. Finally, model 1 can be rewritten in the form
{
N˙ 1 ) -κ(kp11φ1* + kp21(1 - φ1*))N1x2fkd[I] + q1
N˙ 2 ) -κ(kp22(1 - φ1*) + kp12φ1*)N2x2fkd[I] + q2 N˙ I ) -kdNI + qI N˙ S ) qS (9)
where κ ) kt-1/2 and the control variables are the semibatch feed rates of one or both monomers, (q1 and q2) and the feed rate of the initiator (qI). The initiator dynamics are stable and decoupled from the dynamics of N1 and N2. Moreover, the initial state is known with good precision; thus, NI can be estimated in the following manner
N ˆ I ) -kdN ˆ I + qI Finally, as mentioned above, since the overall termination rate constant kt is difficult to calculate theoretically, we decided to estimate its value during the course of a polymerization. The control model used below thus takes on the form
{
N˙ 1 ) -κ(kp11φ1* + kp21(1 - φ1*))N1x2fkd[I] + q1
N˙ 2 ) -κ(kp22(1 - φ1*) ) kp12φ1*)N2x2fkd[I] + q2 N˙ I ) -kdNI + qI k˙ t ) (t) (10)
where (t) is unknown but has bounded dynamics. Recall that we use the heat of reaction to estimate N1 and N2. Thus, the relationship between the number of moles of monomer and the heat of reaction is reflected in the output y(t), given by
y)
∫QR dt - N1tot∆H1 - N2tot∆H2 )
-N1∆H1 - N2∆H2 (11)
3. Application to Solution Copolymerization In this section, we will show how the techniques of section 2 can be applied to the control of the heat generation rate and the copolymer composition rate. Generally, such controls need the on-line measurement of state components (i.e. N1, N2, NI, and NS). The reactor system that we will consider in this work is shown schematically in Figure 1. It consists of a semibatch stirred tank reactor with a cooling jacket and a temperature-controlled thermostatic bath and two inlet streams with separately controlled feed pumps. It is assumed that the reactor temperature TR (continuously monitored by thermocouple TC1) is maintained constant
4818
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
at all times by means of a temperature controller linked to the flow rate at the temperature of the fluid in the cooling jacket. We will suppose that this control loop functions perfectly and not consider it any further in the current work. While it is clear that this last assumption is not perfectly realistic, we have demonstrated that4,20,21 the calorimetric method is robust under nonisothermal conditions. We can therefore be confident that, even if we impose isothermal conditions in the examples presented here (for the sake of simplicity), the results can be extended to nonisothermal systems with imperfect temperature control. It is further assumed that this controller is totally independent of the control system used to regulate the monomer and initiator flow rates that intervene in the control of the composition and the rate of heat generation. Finally, we assume that sufficient temperature and on-line gravimetric data are available to allow us to obtain the rate of heat generation due to reaction as a function of time using an adaptive calorimetric technique.3,4 In the first subsection, we will give an adaptive highgain estimator permitting us to estimate N1, N2, and the variable κ. In the second one, we synthesize nonlinear control laws that will allow us to control both the heat generation rate QR and the copolymer composition. In practice, these controllers used the estimated states given by the high-gain adaptive estimator. 3.1. High-Gain Adaptive Estimator. Our objective here is to estimate the number of moles of both monomers and the variable κ, which depends on N1, N2, and the variable ω. To achieve this, we consider the following augmented system:
()(
)()
N˙ 1 κx2fkd[I]W1 q1 N˙ 2 ) κ 2fk [I]W + q2 x d 2 κ˘ 0
(12)
where the output is defined in eq 11 above, (t) is an unknown function, and the Wj’s (j ) 1, 2) are given by
( ) (
W -(kp11φ1* + kp21(1 - φ1*))N1 W(N) ) W1 ) -(kp22(1 - φ1*) + kp12φ1*)N2 2
)
(13)
To apply the formalism of the high-gain observer recalled in the previous section,3,25-27,37 we consider the following transformation (which is only a slight modification of φ defined in eq A2):
N ) (N1, N2) f (φ1(N), φ2(N)) ) (z1, z2) ) z (14) where
{
z1 ) φ1(N) ) y ) -∆H1N1 - ∆H2N2 z2 ) φ2(N) ) ∆H1W1(N) + ∆H2W2(N)
(15)
Then the first two equations of system 12 can be rewritten in the new system of coordinates (z1, z2) in the following form:
z3 )
(
) ( ) (
0 0 κx2fkdi[I] z+ + κφ0(z) 0 0 q1 0 0 φ1(z) φ2(z) q2
)( )
(16)
where φ0, φ1, and φ2 are some nonlinear functions which depend on z. This implies that if κ is known and constant, then system 12 (the two first equations) is uniformly observable as long as κ(2fkdi[I])1/2 is greater than zero. In this case, an estimator for this system is thus given by system A7 with a1) κ(2fkdi[I])1/2. However, κ is never constant and is difficult to determine exactly, so it should be estimated. If κ is unknown and time varying, then an estimator of system 12 is given by system A10, which can be written
N ˆ˙ ) κˆ (x2fkd[I])W(N ˆ ) + q1g1 + q2g2 + ∂φ -1 T Γ(t)S-1 ˆ 1 + ∆H2N ˆ 2) (17a) θ C (y(t) + ∆H1N ∂N |Nˆ
( )
2
θ1 κˆ˘ ) (y(t) + ∆H1N ˆ 1 + ∆H2N ˆ 2) φˆ
(17b)
with
( ) ( )|
()
()
[ ]
N 2θ 1 0 T , g2 ) , S-1 , N ) N1 , g1 ) θ C ) 0 1 θ2 2 -1 ∂φ1 ∂φ1 1 0 ∂N ∂N ∂φ -1 1 ) ∂φ 1 Nˆ ∂φ 2 Nˆ , Γ(t) ) 0 ∂N Nˆ 2 2 κˆ x2fkd[I] ∂N1 Nˆ ∂N2 Nˆ
[ ] | |
| |
(
)
and finally
φˆ ) (∆H1W1(N ˆ ) + ∆H2W2(N ˆ ))x2fkd[I] Remark: Note that one of the most important characteristics of this estimator is that it does not require any a priori knowledge of (t). Furthermore, and perhaps just as importantly from an operational point of view, calibration of the gain is simple and does not require the resolution of any differential equations (Riccatti equation). 3.2. Control of the Cumulative Copolymer Composition and the Heat Generation Rate. As discussed above, the control scheme should first of all maintain the copolymer composition at a target value (that can vary as a function of time) and, if possible, maintain the highest possible rate of polymerization. In real terms, this means that we must fix a profile for the overall or cumulative composition, which we will manipulate by controlling the flow rate of one (or both) of the monomers. The second point means that we must control the rate of heat generation as closely as possible to the maximum heat removal capacity of the reactor, which will be done by manipulating the initiator flow rate to control the total radical concentration in the reactor. These schemes will be tested under simulation by imposing a profile for QR (usually decreases slightly as a function of time because the overall heat-transfer coefficient will decrease as the quantity of polymer, and thus the viscosity, increases) that represents the input from an adaptive calorimetric measurement scheme.4 To achieve the goals we have set for ourselves, we will use the partial linearization that we reviewed in subsection 3.2. Initially, we will apply this formulation to the tracking of the target value of the copolymer composition using system 9 and the feed rates of the more reactive monomer (all of the less reactive monomer
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4819
is put in the reactor at time t ) 0). The second objective is to optimize the productivity using the feed rates of the initiator. For both cases we need to estimate N1, N2, and kt. It is important to reiterate here that calculating QR from eq 6 is admittedly a simplification. However, as we pointed out above, it has been shown that the calorimetric method can be applied even in large pilot scale reactors. We can therefore be sure that under more realistic conditions we can obtain a good estimate of QR and that for this reason it is acceptable to estimate it from the simplified expression (eq 6). The control model that we will use to synthesize our controller takes the following form
{
x˘ ) g0(x) + u1g1(x) + u2g2(x) H1(x) F ) Ni z˘ ) H(x) ) H2(x) I
( ) ( )
where
()
N1 x ) N2 , NI 0
(
-κ(kp11φ1* + kp21(1 - φ1*))N1x2fkd[I]
g (x) ) -κ(kp22(1 - φ1*) + kp12φ1*)N2x2fkd[I] -kdNI 1 g1(x) ) 0 , g2(x) ) 0
()
)
,
() 0 0 1
u1 ) q1, u2 ) qI, and Fi is the instantaneous composition defined above. In our case the two relative degrees are both equal to 1, and the general form of the control law in eq A13 takes the form
q1 )
1 [-Lf(H1(x)) + v1(t)] Lg1(H1(x))
(18)
where Lg1(h1(x)) and Lf(h1(x)) are nonlinear functions of N1, N2, and kt in which these “variables” are replaced by their estimated values. To track a predefined trajectory of the instantaneous copolymer composition Fd1(t), the control v1(t) is given by
ˆ 1) + v1(t) ) kc1(Fd1 - F
∫
1 (Fd1 - F ˆ 1) dt τ
(19)
where Fd1 is the instantaneous composition set point (can vary as a function of time) for monomer 1 and kc1 and τ are two control tuning parameters. To control the rate of heat generation, we specify the maximum rate of heat generation acceptable (can be a function of time) and use eqs 3 and 6 to identify the desired number of moles of initiator NdI. Our control law then takes the form
qI )
1 [-Lf(H2(x)) + v2(t)] Lg2(H2(x))
(20)
It is easy to verify that in this case we have a linear control law in which Lg2(h2(x)) ) 1 and Lf(h2(x)) ) -kdNI.
And similarly to the first control law, v2(t) is given by
v2(t) ) kr2(NdI - NI) +
∫
1 (NdI - NI) dt τ2
(21)
where kr2 and τ2 are two control tuning parameters. The controllers’ constants are chosen to maximize the speed of convergence to the set point but also to avoid saturating the actuators (here, the pumps). These parameters were set at kc1 ) 0.01, τ ) 105 s, kr2 ) 0.005, and τ2 ) 2 × 108 s. 4. Results and Discussion State Estimation and Control of the Polymerization Reactions. In the current work, the use of the nonlinear observer defined above is demonstrated, and the control formalism is tested using simulated data. To do so, the systems of eqs 6 and 9 were used to simulate the copolymerization of butyl acrylate (BuAs monomer 1) and vinyl acetate (VAcsmonomer 2) in an ethyl acetate solvent (EAc) with benzoyl peroxide (BPO) as the free radical initiator. The rate constants are those shown in Table 1. The initial charge to the reactor was taken as N1,0 ) 0.7 mol, N2,0 ) 6.2 mol, NI0 ) 0.1 mol, and NS,0 ) 34.1 mol. The reaction was assumed to take place under isothermal conditions at 333 K (60 °C). A constant target value for the overall composition of the copolymer was set at a cumulative mole fraction of BuA in the copolymer equal to 0.5. We have chosen here to keep the value of copolymer composition constant, but this is by no means necessary. Note that these initial conditions do not satisfy the composition criterion of FBuA ) 0.5. We chose to start from an erroneous set of initial conditions to show that the observer/control scheme is capable of overcoming such process errors. As stated above, the observer described here uses the instantaneous heat of reaction as an input. To simulate real experiments, QR was calculated from eq 6 using the real values of the charge to the reactor. These real values are those specified above (at time zero) plus the value of any BuA added during the course of the simulated experiment. We imposed a Gaussian noise of zero mean and a standard deviation of (5% on the resulting value of QR. This noisy value is used as the input to the estimator, which then calculates N1, N2, and kt. These estimated values are then used in the control laws. The estimator was provided with initial values of 1010 L‚s/mol and 0.75 and 6.4 mol for kt and the initial number of moles of BuA and VAc in the reactor, respectively. The gains were set at θ ) 0.000 01 and θ1 ) 0.004. Example 1. Composition Control Only. In the first series of simulations, the output from the observer was used to calculate the control of the flow rate of BuA (monomer 1) necessary to maintain the specified value of copolymer composition. No effort was made to control the heat of polymerization. The value of the copolymer composition found using the control law in eq 28 is compared to the set point value in Figure 2. The BuA flow rate needed to achieve this control is also shown in Figure 2. Note that the flow rate of BuA is zero during the initial instants of the reaction, since the conditions are such that we are above the set point composition and there is too much BuA. The rate of heat generation calculated from the conditions in the reactor is shown in Figure 3 for this case. This is the calculated output using the model, not the noisy input to the estimator.
4820
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
Table 1. Data Used in Kinetic and Physical Models property (kg/m3)
density homopolymer density (kg/m3) kp (dm3/(mol s)) kp/kt0.5 ∆Hp (kJ/mol) r1 kd (s-1) density of EAc (kg/m3)
vinyl acetate
butyl acrylate
FVAc ) 1325.3 - 1.34T FBuA ) 1203.8 - 1.04T (K)29 FPVAc ) 121929 FPBuA ) 1372.4 - 1.1T (K)29 independent of temperature 2.7 × 108 exp(-3347/T) T [)] K30 1.58 × 107 exp(-2080/T) T [)] K31 0.40[BuA + VAc] + 0.2632 0.01[BuA + Vac] + 0.1133 89.134 7734 0.02635 5.935 kd ) 4.403 × 1013 exp(-14750/T)32 FVAc ) 1272.6 - 1.27T (K)29
Figure 2. Modeled/estimated (s) copolymer compositions and set point value (- - -), and flow rate of BuA to maintain constant composition.
Figure 3. Rate of heat generation (instantaneous) in semibatch example 1.
The total amount of polymer produced using this technique was approximately 498 g at the end of 4 h (240 min). Example 2. Simultaneous Composition and Rate Control. In the second series of simulations, we chose to simultaneously control the copolymer composition at the same value as specified above, as well as the heat generation rate starting with the same set of initial conditions as in the first example. As discussed above, this allows us to maintain a higher productivity during the course of the semibatch run. To do so, we use the initiator feed flow rate (initiator is dissolved in a solution of EAc at 10% w/w) to alter the rate of reaction. Since the viscosity of the solution in the reactor will increase as the reaction progresses, it is normal for the heat-transfer coefficient of the reactor to decrease slightly as a function of time. Note that although we do add solvent here, it enters the reactor in the initiator stream in such small quantities that it has no effect on heat transfer. Also, although the addition of monomer will of course dilute the reaction system to a certain
(K)29
Figure 4. Heat generation rates and initiator flow rates for example 2: simultaneous semibatch control.
extent, it is added at a rate lower than that of polymer production. For this reason the viscosity will increase as a function of time, thereby lowering the overall heattransfer coefficient. Therefore, without an exact model for solution density, it is not possible to say whether the heat removal capacity of the reactor will increase or decrease with time. We arbitrarily imposed the heat generation rate profile shown in Figure 4 but could just as easily have imposed an increasing profile should we have wished to. In any event, the shape of the maximum heat removal profile is not important here, since we only wish to demonstrate a method of reactor control. Using the same noise level as above and the same model equations, we obtain the “observed” heat generation profile shown in Figure 4. Also shown in Figure 4 is the initiator flow rate needed to maintain the heat profile. Note that we obtain what are effectively occasional shots of initiator added to the reactor instead of a steady, diminishing flow because of the noise in QR and because we are working with a zero-hold operating at 150 s intervals. This allows us to work under more realistic conditions, avoiding situations where we would continuously added “microgram” shots every few seconds (which is not practical, and even if it were, it would be very hard on the pumps). The first addition of initiator in solution does not begin until approximately 300 s. This is not because we are at the maximum heat generation rate during this period but rather because we wish to be sure that the state estimator converges before we begin adding the initiator. Furthermore, we assumed that the initiator feed pump operated at a maximal value of 9 cm3/s. This can also be seen in Figure 3, where the initiator flow rate “flattens out” at 5 × 10-5 mol/s. While we have not discussed the issue here, it is important not to add too much initiator, as this might have a detrimental effect on the molecular weights of the polymer for instance. It remains to be seen what the ideal maximum initiator flow rate should be. The BuA flow rate and copolymer composition are identical to those shown in Figure 2.
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4821
Figure 5. Cumulative quantities of polymer produced with no control over heat of reaction (example 1) and when the heat of reaction was maintained at a safe maximum limit (example 2).
with precision throughout the course of a reaction, especially if one does not have a lot of kinetic data. Also, in cases where the gel or Tromsdorff effect causes the rate constant to vary as a function of time, or if there are significant solvent effects, the fact that we can accurately estimate the value of the termination rate constant as a function of the progress of the reaction without having to rely on a model that might be very inaccurate is a big plus. Thus, the rapid convergence and proven precision of the estimate found here provides us with information that is quite valuable in calculating molecular weights and so forth. Although we have not spoken of the molecular weight distribution, it is important to point out that this quantity can be quite sensitive to the instantaneous free radical concentration and thus the amount of initiator we add. If this turned out to pose a quality problem, the rate of initiation would be constrained by the molecular weight profiles we require, thereby adding a further dimension of difficulty to this problem. This will be treated in a future case study. 5. Conclusions
Figure 6. Variation of the overall termination rate constant with time.
In conclusion, after reviewing the principals of nonlinear observers and control, we looked at their application to a nonlinear system of some economic importance: that is, a free radical copolymerization. The results discussed in this current paper show how to develop control laws for systems such as those discussed elsewhere.4 We have demonstrated their feasibility and the fact that such nonlinear systems both converge rapidly and are relatively robust to poor inputs and process noise. Nomenclature
As stated above, the composition problem is uncoupled from the heat generation problem. It can be seen in Figure 4 that the total amount of heat removed from the reactor is greater in the second example than the first. In fact, maintaining the highest possible rate of heat generation allows us to produce approximately 567 g of polymer compared to only 498 g for the first example. A comparison of the total amount of polymer produced in both cases is shown in Figure 5. The fact that we control the rate of heat generation allows us to produce about 12% more polymer. Note that the heat generation rate profile imposed here is rather conservative. It is possible, at least when using simulated data, to impose any profile one wishes. It can be seen that there is satisfactory agreement between estimated and simulated data in Figure 6. Also, despite the very high value for kt supplied to the observer at time zero, it can be seen that the estimated value of the overall termination rate constant shown in Figure 6 also converges toward the theoretical value calculated from estimator input concentrations and a previously presented model.28 Note that the speed of convergence would be improved if we had supplied the estimator with an initial value of the overall rate constant that was much closer to the “true” value. As can be seen in Figure 6, the initial guess for kt was in fact 2 orders of magnitude greater than the “real” value. In theory it should be possible to obtain a much better estimate of kt0 for most systems by simple experimentation, but we chose to demonstrate the robustness of the estimator to very poor initial conditions. Nevertheless, even if it is possible for us to estimate kt to within an order of magnitude, it remains difficult to calculate it
cpr ) heat capacity of reactor contents and fittings f ) initiator efficiency fi ) mole fraction of species i relative to total number of moles of all monomers Fi ) mole fraction of i in polymer ∆Hp,i ) heat of polymerization of species i [I] ) concentration of initiator kd ) initiator decomposition rate constant kpij ) propagation rate constant of radical of type i with monomer of type j kt ) total termination rate constant mr ) mass of reactor contents, fittings, etc. Nx, x ) 1, 2, S, I ) number of moles of monomers 1 and 2 and of solvent and of initiator Nxtot ) total number of moles of species x fed to reactor qx,in, x ) 1, 2, S, I ) molar flow rate of species x into reactor r1, r2 ) reactivity ratios of monomers 1 and 2 [R*] ) total concentration of free radicals in reactor Rpi ) rate of polymerization of monomer i Vp ) volume of reactor contents xi ) individual conversion of monomer i Greek Symbols φi* ) mole fraction of radicals terminating in radical of type i ω ) adjustable parameter for calculating termination rate coefficient ()1 here)
Appendix. A Short Review of Nonlinear Observers and Geometric Control A.1. Nonlinear Observers. In this subsection, we restrict ourselves to the class of single output systems:
4822
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
{
y ) h(x)
LX(φ)(x) ) φ(x + X(x)) - φ(x) lim
m
)
f0
∂φ
Xj(x) (x) ∑ ∂x j)1 j
LkX(φ),
Xj is the jth component of X. By we denote the kth Lie derivative of φ with respect to X, LkX(φ) ) 0 LX(Lk-1 X (φ)) with the following convention: LX(φ) ) φ. Now, set T (A2) z ) φ(x) ) [h(x), Lg0(h)(x), ..., Lgk-1 0 (h)(x)]
and assume that φ is one to one and ∂φ/∂x(x) is nowhere zero. It has been shown29,30 that if system A1 is uniformly observable, then φ transforms A1 into the following canonical form:
y ) Cz
ujψ (z) ∑ j)1
( ) 0 1 . . .
with
1 (s,u,z1, ..., zn-1), ψ1n(s,u,z)]T ψn-1
Under a similar hypothesis, as in ref 35, if, moreover,
(H1) 0 < F1 e ai(s(t)) e F2, ∀t g 0 for some constants F1, F2
xˆ3 ) g0(xˆ ) +
ujgj(xˆ ) ∑ j)1
() ∂x
-1
|
d (a (s(t))) e M dt i
(
)
where Sθ is as above and
1 0 L) . 0 0
for j ) 1, ..., m
T S-1 ˆ - y) j C (Cx
|
then an exponential observer for eq A5 takes the following form
(
This canonical form is then used to construct an exponential observer for system A1 under some additional hypotheses.37 It is shown that an exponential observer for system A1 takes the form
∂φ
)
a1(s) 0 . 0 a2(s) . 0 0 , C ) (1, 0, ..., 0), . . . 0 . . . an-1(s) . . . 0 0 ψ (s,u,z) ) [0, ..., 0, ψ0n(s,u,z)]T, and ψ1(z) ) [ψ11(s,u,z1), ψ12(s,u,z1,z2), ...,
j (z1, ..., zn-1), ψj(z) ) [ψj1(z1), ψj2(z1,z2), ..., ψn-1
m
(A6)
-1 ∂φ T xˆ3 ) g(s,u,xˆ ) (xˆ ) Λ-1S-1 ˆ ) - y) θ C (h(x ∂x
. 0 . 0 . 0 , C ) (1, 0, ..., 0), . 1 . 0 ψ0(z) ) [0, ..., 0, ψ0n(z)]T, and ψjn(z)]T
(
0 0 A) . 0 0
(H2) ∃M > 0; ∀t g 0,
where
1 0 . . .
{
z3 (t) ) A(s(t)) z(t) + ψ0(s(t), u(t),z(t)) + ψ1(s(t), u(t), z(t)) y(t) ) Cz(t)
j
(A3)
(A5)
where x(t) ∈ Rn, u(t) ∈ Rm, y(t) ∈ R, and s(t) is a smooth signal from R+ into Rk. Assume that there exists a change of coordinates z ) φ(x) which transforms A5 into the following form
m
0
z3 ) Az + ψ (z) +
0 0 A) . 0 0
x3 (t) ) g(s(t), u(t), x(t)) y(t) ) h(x)
(A1)
x(t) ∈ Rn, ui(t) ∈ R, u ) (u1, ..., um), and y(t) ∈ R. System A1 is said to be uniformly observable if, for every two initial states x * x and every two admissible inputs u defined on any [0, T], there exists t ∈ [0, T] such that y(x,u,t) * y(x,u,t), where y(x,u,t) is the output associated to the initial state x and the control u. Now, consider any differentiable function from Rn into R and set X any vector field on Rn (i.e. a differentiable map from Rn into Rn). The Lie derivative of φ with respect to X denoted by LX(φ) is defined by
{
{
m
uigi(x) ∑ i)1
x3 ) g0(x) +
(A4)
0 a1(s) . . .
0 0 a1(s) a2(s) . .
. . . . 0
0 0 0 0 a1(s) ... an-1(s)
(A7)
)
It has been demonstrated in the literature38 that eq A7 converges. Consider now the following nonlinear system:
{
x3 (t) ) µg0(s(t), x(t)) + g1(s(t), u(t), x(t)) y(t) ) h(x)
(A8)
θSθ + ATSθ + SθA ) CTC
where x(t) ∈ Rn, u(t) ∈ Rm, y(t) ∈ R, and s(t) is a smooth signal from R+ into Rk. Assume that there exists a change of coordinates z ) φ(x) which transforms eq A8 into the following form:
where θ is a positive constant parameter. In what follows, we well extend this observer synthesis to more general systems. Consider a nonlinear system:
z3 (t) ) µAz(t) + ψ0(s(t), u(t), z(t)) + ψ1(s(t), u(t), z(t)) y(t) ) Cz(t)
where Sθ is the unique solution of
{
(A9)
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4823
with
( )
0 0 A) . 0 0
φn(x)) becomes a local system of coordinates and system A8 takes the following normal form:
1 0 . 0 0 1 . 0 . . . 0 , C ) (1, 0, ..., 0), . . . 1 . . . 0 ψ0(s,u,z) ) [0, ..., 0, ψ0n(s,u,z)]T, and ψ1(z) ) [ψ11(s,u,z1), ψ12(s,u,z1,z2), ..., 1 (s,u,z1,...,zn-1), ψ1n(s,u,z)]T ψn-1
In the case where the parameter µ is known and constant, system A9 is a particular case of system A6. This implies that system A7 is an exponential observer of eq A9 in which a1 ) a2 ) ..... ) an-1 ) µ. However, if the parameter µ is unknown and time varying, then we propose the following nonlinear adaptive observer:
xˆ3 ) µg0(s,u,xˆ ) + g1(s,u,xˆ ) -1 ∂φ T (xˆ ) Λ-1S-1 ˆ ) - y) θ C (h(x ∂x
(
)
Set
2
θ (h(xˆ ) - y) µˆ˘ ) φ(xˆ )
(A10)
where φ(xˆ ) is a nonlinear function which depends on x. A.2. Nonlinear Geometric Control. In this subsection, we will review some well-known general results concerning the input-output feedback partial linearization. More details can be found in the literature.39,40 Consider the multiinput, multioutput, nonlinear system
{
m
x3 ) g0(x) +
where
()
uigi(x) ∑ i)1
(A11)
y ) h(x)
()
x1 u1 ∈ Rm, and x ) : ∈ R n, u ) : xn um
Ω(z) )
and
(
( )
Lgr10(h1(z)) b(z) ) Lgr20(h1(z)) Lgrm0 (h1(z))
Lg1Lgr10-1(h1(z)) .... LgmLgr10-1(h1(z))
(
Lg1Lgrm0 -1(hm(z)) .... LgmLgrm0 -1(hm(z))
1 k11 ... kr1 0 ...
k21
K ) 0 ... 0 0 ... 0
...
... 0 kr22
0 ... km 1
0 ... ...
0 m ... krm
)
)
Then, the following control law
( )
h1(x) ∈ Rm h(x) ) : hm(x)
We suppose for this system that the distribution ∆ spanned by ∆ ) {g1 .... gm} is involutive (which means that the Lie bracket [gi gj] ) (∂gj/∂x)gi - (∂gi/∂x)gj ) m Rl(x)gl). This hypothesis is easily verified for our ∑l)1 system, since the vectors gl, ..., gm are constant. Recall that the relative degrees {rl, ..., rm} of eq A8 are the smallest integers ri g 1 such that for all x in a neighborhood of x0, for all j, 1 e j e m, and for all i, 1 e i e m, we have LgjLgri0-1(hi(x0)) * 0 for at least one j, 1 e j e m. Set r ) rl + ... + rm and
φr ) [h1, Lg0(h1), ..., Lgrl0-1(h1), ..., hm, Lg0(hm), ..., Lgrm0 -1(hm)]T ) [φ1, ..., φr]T and assume that (∂φr/∂x)|x)x0 is of rank r. If r is strictly less than n, it is always possible to choose n - r functions φr+1, ..., φn such that x f z ) φ(x) ) (φ1(x), ...,
u(t,z) ) Ω-1(z)[-b(z) - KOr + v(t)]
(A13)
gives rise to the following open-loop system:
{
z3 1 ) Az1 + Bv(t) z3 2 ) ψ(z) y ) Cz
(A14)
1
where v is the new control variable that is used for the tracking of the desired trajectory. A is a diagonal block matrix given by
(
A1 0 A) : 0
0 A2 : 0
... ... : ...
0 0 : Am
)
B ) (er1, ..., erm), and C ) (e1, er1+1, ..., erm-1+1), where ei are vectors of dimension r in which all the components vanish except the ith component, which is equal to one. The elements of block A are
4824
( )
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
0 0 Ai ) . 0 ki1
1 0 . . .
0 1 . . .
. . . . .
( )
0 0 v1(t) 0 , v(t) ) : , z1 ) [z1, ..., zr]T 1 vm(t) kiri
and z2 ) [zr+1, ..., zn]T. The gains ki1 ... kri i are chosen such that the real parts of the eigenvalues of Ai are strictly negative. The local asymptotic stability of eq A9 is conditioned by the stability of the zero dynamics: z3 2 ) ψ(0,z2).39,40 Literature Cited (1) Penlidis, A.; Ponnuswamy, S.; Kiparissides, C.; O’Driscoll, K. Polymer reaction engineering: modelling considerations for control studies. Chem. Eng. J. 1992, 50, 95. (2) McKenna, T.; Fe´votte, G.; Graillat, C.; Guillot, J. Joint Use of Calorimetry, Densimetry and Mathematical Modelling for Multiple Component Polymerisations. Chem. Eng. Res. Des. 1996, 74A, 340. (3) Fe´votte, G.; McKenna, T.; Othman, S.; Hammouri, H. Nonlinear tracking of glass transition temperatures for free radical emulsion copolymers. Chem. Eng. Sci. 1998, 53, 773. (4) Othman, S.; Barudio, I.; Fe´votte, G.; McKenna, T. On-line monitoring and modelling of free radical copolymerisations: Butyl Acrylate/Vinyl Acetate. Polym. React. Eng. 1999, 7 (1), 1. (5) Kiparissides, C.; Morris, J. Intelligent manufacturing of polymers. Comput. Chem. Eng. 1998, 20, S1113. (6) Dimitratos, J.; Elic¸ abe, G.; Georgakis, C. Control of Emulsion Polymerization Reactors. AIChE J. 1994, 40, 1993. (7) Congalidis, J. P.; Richards, J. R. Process control of polymerization reactors: An industrial perspective. Polym. React. Eng. 1998, 6, 71. (8) Kozub, D. J.; MacGregor, J. F. Feedback control of polymer quality in semi-batch copolymerization reactors. Chem. Eng. Sci. 1992, 47, 1047. (9) Mutha, R. K.; Cluett, W. R.; Penlidis, A. A new multirate measurement-based estimator: Emulsion copolymerization batch reactor case study. Ind. Eng. Chem. Res. 1997, 36, 1036. (10) Mutha, R. K.; Cluett, W. R.; Penlidis, A. Nonlinear modelbased estimation and control of a polymer reactor: Experimental results. AIChE J. 1997, 43, 3042. (11) Dochain, D.; Pauss, A. On-line estimation of microbial specific growth rates: an illustrative case study. Can. J. Chem. Eng. 1988, 47, 327. (12) Kravaris, C.; Wright, R. A.; Carrier, J. F. Nonlinear controllers for trajectory tracking in batch processes. Comput. Chem. Eng. 1989, 13, 73. (13) 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. (14) Schmidt, C. U.; Reichert, K. H. Reaction Calorimeter. A contribution to Safe Operation of Exothermic Polymerizations. Chem. Eng. Sci. 1988, 43, 2133. (15) Moritz, H. U. Polymerization Calorimetry: A powerful Task for Reactor Control. In 3rd International Workshop on Polymer Engineering; Reichert, K. H., Geisler, U., Eds.; VCH Verlag: Berlin, 1989. (16) Urretabizkaia, A.; Leiza, J. R.; Asua, J. M.; On-line terpolymerization control by means of a semi-continuous emulsion copolymerization. AIChE J. 1994, 40, 1850. (17) Saenz de Buruaga, I.; Armitage, P. D.; Leiza, J. R.; Asua, J. M. On-line Control for Maximum Production Rate of Emulsion Polymeris of Well-Defined Polymer Composition. In Proceedings of the First European Congress on Chemical Reaction Engineering, Florence, Italy, 1997; Vol. 1, p 117. (18) Saenz de Buruaga, I.; Arotc¸ arena, M.; Armitage, P. D.; Gugliota, L. M.; Leiza, J. R.; Asua, J. M. On-line calorimetric control of emulsion polymerization reactors. Chem. Eng. Sci. 1996, 51, 2781. (19) Saenz de Buruaga, I.; Armitage, P. D.; Asua, J. M. Nonlinear Control for Maximum Production Rate of Latexes of Well-Defined Polymer Composition, Ind. Eng. Chem. Res. 1997, 36, 4243.
(20) Fe´votte, G.; McKenna, T. F.; Othman, S.; Santos, A. M. A combined Hardware/software sensing approach for on-line control of emulsion polymerisation processes. Comput. Chem. Eng. 1998, 22, S443. (21) McKenna, T. F.; Othman, S.; Fe´votte, G.; Santos, A. M. Integrated approach to monitoring, state estimation and control of polymer reactors. In 6th International Workshop on Polymer Reaction Engineering; Reichert, K. H., Moritz, H. U., Eds.; WileyVCH: Berlin, 1998. (22) Gloor, P. E.; Warner, R. J. Developing feed policies to maximize productivity in emulsion polymerization processes. Thermochim. Acta 1996, 289, 243. (23) Hamielec, A.; MacGregor, J.; Penlidis, A. Multicomponent free-radical polymerisation in batch, semi-batch and continuous reactors Makromol. Chem. Macromol. Symp. 1987, 10/11, 521. (24) Dube´, M. A.; Soares, J. B. P.; Penlidis, A.; Hamielec, A. E. Mathematical Modelling of Multicomponent Chain-Growth Polymerizations in Batch, Semi-batch and Continuous Reactors: A review. Ind. Eng. Chem. Res. 1997, 36, 966. (25) Deza, F.; Busvelle, E.; Gauthier, J. P.; Rakotopara, D. High gain estimation for nonlinear systems. Syst. Control Lett. 1992, 18, 295. (26) Gibon-Fargeot, A. M.; Hammouri, H.; Celle, F. Nonlinear estimators for chemical reactors. Chem. Eng. Sci. 1994, 49, 2287. (27) Farza, M.; Hammouri, H.; Othman, S.; Busawon, K. Nonlinear observers for parameter estimation in bioprocesses. Chem. Eng. Sci. 1997, 52/53, 4251. (28) McKenna, T. F.; Fortuny, M. Effect of solvent on the rate constants in the solution copolymerisation of butyl acrylate and vinyl acetate. Polym. J. 1999, 31, 7. (29) Barudio, I.; Fe´votte, G.; McKenna, T. F. Density data for copolymer systems: Butyl Acrylate/Vinyl Acetate Homo- and Copolymerisation in Ethyl Acetate. Eur. Polym. J. (to appear). (30) Hutchinson, R. A.; Richards, J. R.; Aronson, M. T. Determination of propagating rate coefficients by pulsed-laser polymerisation for systems with rapid chain growth: vinyl acetate. Macromolecules 1994, 27, 4530. (31) Lyons, R. A.; Hutovic, J.; Piton, M. C.; Christie, D. I.; Clay, P. A.; Manders, B. G.; Kable, S. H.; Gilbert, R. G. Pulsed-laser polymerisation measurements of the propagation rate coefficient for butyl acrylate. Macromolecules 1996, 29, 1918. (32) McKenna, T. F.; Villanueva, A.; Santos, A. M. Effect of solvent on the rate constants in solution polymerisation. Part I: butyl acrylate. J. Polym. Sci., Part A: Polym. Chem. (to appear). (33) McKenna, T. F.; Villanueva, A. Effect of solvent on the rate constants in solution polymerisation. Part II: vinyl acetate. J. Polym. Sci., Part A: Polym. Chem. (to appear). (34) Busfield, W. K. Heats and Entropies of Polymerization, Ceiling Temperatures, Equilibrium Monomer Concentrations, and Polymerizability of Heterocyclic Compounds. In Polymer Handbook; Brandrup, J., Immergut, E. I., Eds.; Wiley: New York, 1989. (35) Dube´, M. A.; Penlidis, A. A systematic approach to the study of multicomponent polymerization kinetics-the butyl acrylate/methyl methacrylate/vinyl acetate example: 1. Bulk copolymerization. Polymer 1995, 36, 587. (36) Gauthier, J. P.; Bornard, G. Observability for any u(t) of a class of bilinear systems. IEEE Trans. Autom. Control 1981, 26, 922. (37) Gauthier, J. P.; Hammouri, H.; Othman, S. A simple observer for nonlinearsystems, applications to bioreactors. IEEE Trans. Autom. Control 1992, 37, 875. (38) Farza M.; Busawon, K.; Hammouri, H. Simple nonlinear observer for on line estimation of kinetics rates in bioreactors. Automatica 1998, 34, 301. (39) Nijmeijer, H.; Van der Schaft, A. Non linear dynamical control systems; Springer-Verlag: Berlin, 1990. (40) Isidori, A. Nonlinear control systems; Springer-Verlag: Berlin, 1989.
Received for review November 5, 1998 Revised manuscript received August 13, 1999 Accepted August 17, 1999 IE9806996
