A New Multirate-Measurement-Based Estimator: Emulsion

This paper presents the development of a state estimator for systems with multirate measurements. A fixed-lag smoothing based extended Kalman filter ...
0 downloads 0 Views 231KB Size
Download PDF
1036

Ind. Eng. Chem. Res. 1997, 36, 1036-1047

A New Multirate-Measurement-Based Estimator: Emulsion Copolymerization Batch Reactor Case Study R. K. Mutha and W. R. Cluett* Department of Chemical Engineering, University of Toronto, Toronto, Ontario, Canada M5S 3E5

A. Penlidis Department of Chemical Engineering, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1

This paper presents the development of a state estimator for systems with multirate measurements. A fixed-lag smoothing based extended Kalman filter algorithm is proposed. The performance of the algorithm has been evaluated using a case study on a complex emulsion copolymerization batch process. The proposed estimator uses measurements multiple times to achieve good convergence properties and robustness to state and measurement noise. The simulations also demonstrate that the performance of the proposed algorithm is superior to the basic extended Kalman filter. 1. Introduction One of the major obstacles in the control of polymer reactors is the lack of adequate and frequent measurements. In the literature, several estimation schemes based on the extended Kalman filter (EKF) have been proposed for process state estimation and subsequent inclusion in process control schemes in order to attempt to compensate for this lack of measurements. MacGregor et al. (1984) and Elicabe and Meira (1988) provided surveys on the estimation and control of polymerization reactors. Chien and Penlidis (1990) reviewed the impact of slow measurements on polymerization reactor control and the need for using online state estimation. Schuler and Schmidt (1993) have recently summarized the basics of state estimation and its application to polymerization reactors, bioreactors, and other chemical reactors. Dimitratos et al. (1994) have also recently reviewed various estimation and control applications to emulsion polymerization reactors. Kozub and MacGregor (1992) developed a fixed-point reiterated EKF for faster convergence of the states of semibatch polymerization reactors. They also demonstrated the importance of including nonstationary states in the estimator to handle bias and model-process mismatch. Dimitratos et al. (1989) simulated an adaptive hybrid EKF for emulsion copolymerization state estimation. Subsequently, Dimitratos et al. (1991) applied this EKF to an experimental emulsion copolymerization system. Adebekun and Schork (1989) also investigated the performance of an EKF applied to a solution polymerization process and implemented the estimator during feedback control of the reactor at the simulation level. Systems with multirate measurements can be divided into two categories: one with primary and secondary measurements and the second with fast (frequent) and slow (infrequent) measurements. In the first type, secondary inferential measurements are available more frequently than the primary measurements, and the * Telephone: 416-978-6697. Fax: 416-978-8605. E-mail: [email protected]. S0888-5885(96)00100-5 CCC: $14.00

estimator states can be made observable by using secondary inferential measurements even in the absence of primary measurements (Gudi et al., 1995). In the second type of systems, the structure of the process model and, hence, the flow (direction of propagation) of information is such that complete state observability cannot be met without the slow measurements. Polymerization processes usually belong to this second category of multirate systems. For instance, Schuler and Suzhen (1985) have shown that no information flow exists from the different moments of the chain length distribution (quality part of the model) to conversion (productivity part of the model). Therefore, observability for the moments of the chain length distribution cannot be met by using density (conversion) measurements alone. One other point of clarification worth mentioning here is that fast measurements are assumed to be instantaneous in the sense that the measurement is taken and is available at the same sampling instant. On the other hand, slow measurements are taken at one sampling instant but are not available until some later sampling instant due to time delay associated with sample analysis, i.e., gel permeation chromatography (GPC). Very few of the estimation algorithms in the literature are designed to handle systems with multirate measurements. One exception is the work by Ellis et al. (1988) where an EKF was designed for and applied to a batch polymerization system with fast and slow measurements. In their algorithm, the calculated state estimates from the faster measurements are discarded when the slower measurements become available and the estimator is applied from the time the slower measurements are taken. This discarding of the estimated states and restarting of the estimator from the time the slower measurements are taken is inefficient and can lead to slower convergence, since useful information is unnecessarily being ignored. Also, this approach does not address systems having several slow measurements each with a different delay. In this paper, a new extended Kalman filter algorithm with fixed-lag smoothing (as opposed to fixed-point smoothing) for estimating the states of nonlinear sys© 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1037

tems with multirate measurements is presented. The proposed algorithm, which does not discard any useful information and can handle several measurements with different amounts of delay, will be shown to result in better convergence than the basic EKF algorithm with all measurements available without delay. The smoothing of the states can be implemented at all times (including startup), without any increase in computational load over time, which would be the case with a fixed-point smoothing algorithm. The algorithm is also applicable to systems with variable measurement delays. To illustrate the use of the proposed estimator, a case study on an emulsion acrylonitrile/butadiene copolymerization batch system is presented. The organization of the paper is as follows: section 2 briefly summarizes the Kalman filter with fixed-lag smoothing, section 3 presents the new extended Kalman filter algorithm with fixed-lag smoothing for multirate systems, section 4 describes the model of the acrylonitrile/butadiene emulsion copolymerization batch process, and finally, section 5 evaluates the proposed algorithm through numerous simulations on this batch process. 2. Kalman Filtering and Smoothing Smoothing is defined as the estimation of past states based on past but more recent measurements (X ˆ (k i|k) ∀ i g 1, where k is the current sampling instant). Filtering is defined as the estimation of current states based on past and current measurements (X ˆ (k|k)), and prediction is defined as the estimation of future states based on past and current measurements (X ˆ (k + i|k) ∀ i g 1). The main motivation for using smoothing is that smoothed states actually use more data (by using measurements beyond the time of the states being estimated) and, therefore, have better convergence properties than filtered states. However, the smoothing algorithms are related to the filtering problem and are, in fact, based on modifications to the standard Kalman filtering algorithm. An early survey paper on smoothing by Meditch (1973) describes its application to linear and nonlinear dynamic systems. State smoothing can be implemented in three different ways (Anderson and Moore, 1979; Elbert, 1984), each with a different objective. These three smoothing algorithms are fixed point, fixed lag, and fixed interval. All three types of smoothing use the entire sequence of available measurements. Fixed-point smoothing estimates the state vector at a single fixed point in time, i.e., X ˆ (k0|k), where k0 is a sampling instant fixed in time (k0 < k). Fixed-lag smoothing estimates the state vector of the system at a specified time lag from the current sampling instant, i.e., X ˆ (k - SI|k), where SI denotes the smoothing interval. Fixedinterval smoothing estimates the state vector at all sampling instants including the current sampling instant, i.e., X ˆ (i|k) ∀ 0 e i e k. The algorithm proposed in this paper is based on fixed-lag smoothing. Use of fixed-point and fixedinterval smoothing would result in an excessive computational load. In addition, it is well-known that smoothing algorithms in general have a negligible effect on the smoothed states beyond two to three times the dominant time constant of the Kalman filter (Anderson and Moore, 1979). The fixed-lag smoothing algorithm is well documented in the literature, e.g., Anderson and Moore (1979) and Elbert (1984). Therefore, only the key

features of this algorithm are summarized in the next subsection. 2.1. Fixed-Lag Smoothing Algorithm. Consider a time-varying, discrete-time linear dynamic system with n states, m measurements, and p inputs described by the equations

X0(k|k - 1) ) A(k,k - 1)X0(k - 1|k - 1) + B(k,k - 1)U(k - 1) + C(k,k - 1)W(k - 1) Y(k) ) H(k)X0(k) + V(k)

(1)

where W(k) and V(k) are zero mean, white noise sequences with covariances equal to Q(k) and R(k), respectively. The dimensions of the above terms are as follows: A (n × n), B (n × p), C (n × n), H (m × n), Q (n × n), R (m × m), X0 (n × 1), U (p × 1), W (n × 1), Y (m × 1), and V (m × 1). The notation A(k,k - 1), B(k,k 1), and C(k,k - 1) indicates the values of these matrices between sampling instants k - 1 and k. The deterministic U(k) can be ignored without loss of generality. For nonlinear systems, the above set of equations is obtained by linearizing the system at the operating point. Smoothing of the states is achieved by using a state vector augmentation technique, with the number of augmented state vectors being equal to the smoothing interval, SI. The augmented system can be written as

[ ] [ ][ ] [ ] [] X0(k|k - 1) X1(k|k - 1) X2(k|k - 1) ) · · · XSI(k|k - 1) A(k,k - 1) I 0 · · · 0

0 0 I · · · 0

· · · 0 · · · 0 · · · 0 · · · · · · I

0 0 0 · · · 0

X0(k - 1|k - 1) X1(k - 1|k - 1) X2(k - 1|k - 1) + · · X (k -· 1|k - 1) SI

C(k,k - 1) 0 0 W(k - 1) · · · 0

X0(k) X1(k) X (k) Y(k) ) [H(k) 0 0 · · · 0] 2 · + V(k) · · XSI(k)

(2)

The augmented state vectors in eq 2 are related according to

XSI(k) ) XSI-1(k - 1) ) ... ) X1(k - SI + 1) ) X0(k - SI) Therefore, the augmented system contains the SI previous values of the state vector X0(k). Application of the Kalman filter to this augmented system results in smoothed estimates of the states at a fixed lag of SI

1038 Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997

sampling intervals; i.e., XSI(k|k) ) X0(k - SI|k). By doing this, the order of the augmented system has increased SI times. However, the matrices of the augmented system are sparse and result in relatively simple equations. More details of the augmented Kalman filter can be found in Elbert (1984). The two sets of equations used recursively to form the covariance matrices for this smoothing algorithm are given as

estimates with the above Kalman gains are given by the following equations:

P00(k|k) ) [I - K0(k)H(k)]P00(k|k - 1)

X ˆ 0(k - 2|k) ) X ˆ 0(k - 2|k - 1) + K2(k)[Y(k) H(k)X ˆ 0(k|k - 1)]

P10(k|k) ) P10(k|k - 1)[I - HT(k)KT0 (k)]

X ˆ 0(k|k) ) X ˆ 0(k|k - 1) + K0(k)[Y(k) H(k)X ˆ 0(k|k - 1)] X ˆ 0(k - 1|k) ) X ˆ 0(k - 1|k - 1) + K1(k)[Y(k) H(k)X ˆ 0(k|k - 1)]

l (k)KT0 (k)]

T

P20(k|k) ) P20(k|k - 1)[I - H

X ˆ 0(k - SI|k) ) X ˆ 0(k - SI|k - 1) + KSI(k)[Y(k) H(k)X ˆ 0(k|k - 1)] (6)

l PSI,0(k|k) ) PSI,0(k|k - 1)[I - HT(k)KT0 (k)]

(3)

and

P00(k|k - 1) ) A(k,k - 1)P00(k - 1|k - 1)AT(k,k - 1) + C(k,k - 1)Q(k)CT(k,k - 1) T

P10(k|k - 1) ) P00(k - 1|k - 1)A (k,k - 1) P20(k|k - 1) ) P10(k - 1|k - 1)AT(k,k - 1) l PSI,0(k|k - 1) ) PSI-1,0(k - 1|k - 1)AT(k,k - 1)

(4)

where Pi,0 are the covariances of the Xith state(s) with the X0 state(s). The smoothing (Kalman) gain matrices are given by the following set of equations:

K0(k) ) P00(k|k - 1)HT(k)[H(k)P00(k|k - 1)HT(k) + R(k)]-1 K1(k) ) P10(k|k - 1)HT(k)[H(k)P00(k|k - 1)HT(k) + R(k)]-1 K2(k) ) P20(k|k - 1)HT(k)[H(k)P00(k|k - 1)HT(k) + R(k)]-1 l KSI(k) ) PSI,0(k|k - 1)HT(k)[H(k)P00(k|k - 1)HT(k) + R(k)]-1 (5) Note that the first equation in each of the above sets of equations (3), (4), and (5) corresponds to the standard Kalman filtering equations without any smoothing. Equation set (5) indicates that there is only one matrix to be inverted, and it is the same matrix that would have to be inverted even without augmentation. Therefore, the computational load due to matrix inversions has not increased for the augmented system. The smoothed

where X ˆ 0(k|k - 1) is computed using the process model equation (1). The above implementation results in the estimation of the states at all lags less than or equal to SI. There are other approaches to the fixed-lag smoothing problem that yield only state estimates at lag SI. 3. EKF with Fixed-Lag Smoothing for Multirate Systems This section develops the proposed extended Kalman filter algorithm with fixed-lag smoothing for multirate systems, where the term “extended” refers to the fact that the algorithm is applied to a linearized version of the nonlinear system. The proposed algorithm is based on a “boot-strap” application of filtering and smoothing. The value of the lag used for smoothing is a key parameter of this algorithm and is referred to as the smoothing interval (SI). Two EKF’s are used in this scheme: one for filtering the states (EKF1) based on the process model and the second for smoothing the states (EKF2). EKF1 is used for estimating the states at time k with the available measurements at time k, i.e., for calculating X ˆ (k|k). The standard extended Kalman filter is used for EKF1, and the states associated with EKF1 are denoted by X ˆ 1. The first equation from each set of equations (3), (4), (5), and (6) corresponds to the standard EKF algorithm and is used by EKF1. EKF2 is used for smoothing the states, i.e., for calculating X ˆ (k - SI|k). The complete sets (3), (4), (5), and (6) compose the fixed-lag smoothing algorithm and are used by EKF2. The states associated with EKF2 are denoted by X ˆ 2. In general, the states of EKF2 are a subset of the states of EKF1; i.e., it is often unnecessary to smooth all of the states of X ˆ 1. For multirate systems, a different number of measurements is available at consecutive sampling instants. The EKF can easily accommodate this change in size of the measurement vector by using only the appropriate rows of H(k). Corresponding to these rows of H(k), a submatrix of the noise covariance matrix R(k) is used by the estimator. Then, the Kalman gain matrix is computed with the varying (in size) H and R matrices, where the number of columns in the Kalman gain matrix is equal to the number of measurements. The proposed algorithm assumes that all states in both EKF1 and EKF2, with X ˆ 1 and X ˆ 2 considered separately, are observable. Since the algorithm is based on the extended Kalman filter equations, it is an application of a linear algorithm assuming piecewise linearity of the nonlinear system. In general, linear

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1039

Figure 1. Fixed-lag smoothing-based EKF for systems with multirate measurements. yf are the fast measurements, ys are the slow measurements, EKF1 is the filter, and EKF2 is the smoother.

filters implemented on nonlinear systems can diverge, even though observability conditions are met. 3.1. Proposed Algorithm. The convergence of the estimated states, which are observable only with the slow measurements, to their true values depends on how frequently slow measurements are taken and on their respective measurement delay. In polymerization processes, measurements are often associated with long measurement delays and are taken infrequently. The convergence properties of these states can be enhanced by using the available measurements several times. In other words, a slow measurement should be used at each and every sampling interval until the next slow measurement becomes available, along with the new fast measurements, to estimate the states. The proposed algorithm tries to compensate for infrequent slow measurements by repeated use of the available slow measurements. For the sake of clarity in presenting the algorithm, a two-rate measurement system will be considered where the first set of measurements is available instantaneously and the second set is taken at a slower rate and is available with a delay. It is also assumed for simplicity that all of the slower measurements become available with the same amount of delay equal to d sampling intervals. Consider that at sampling instant k, the following variables are available from the previous sampling ˆ 1(k|k - 1), X ˆ 2(k - i|k instants: X ˆ 1(k - SI|k - 1), K1(k), X 2 1 - 1) ∀ 0 e i e SI, P (k|k), and Pi,0(k|k) ∀ 0 e i e SI. The implementation of the proposed algorithm is as follows: 1. At sampling instant k, the fast (instantaneous) measurements (yf(k)) are used for filtering the predicted states of EKF1; i.e.,

X ˆ 1(k|k) ) X ˆ 1(k|k - 1) + K1(k)(yf(k) - yˆ f(k)) 2. This step is executed only when the slow measurements are taken. Calculate the entire Kalman gain matrix (K2i (k) ∀ 0 e i e SI) of EKF2. For the slow measurements, the Kalman gain columns are stored in 2′ a buffer, K2′ i (k). This buffer is referred as Ki (k - d) after d sampling instants when these slow measurements become available. 3. This step is executed only when slow measurements (ys(k - d)) become available. Update the back end of the smoothed states (EKF2) as follows:

X ˆ 2(k - i|k - 1) ) X ˆ 2(k - i|k - 1) + 2′ (k - d)(ys(k - d) - yˆ s(k - d)) ∀ d e i e SI Ki-d

4. Using yf(k), EKF2 smooths the states to calculate X ˆ 2(k - SI|k). This involves the execution of the sets of equations (4), (5), (6), and (3) for X ˆ 2. 5. The smoothed states X ˆ 2(k - SI|k) are used to update the EKF1 states from X ˆ 1(k - SI|k - 1) to X ˆ 1′(k - SI|k). 6. The states X ˆ 1′(k - SI|k) are updated using all measurements available at sampling instant k - SI to generate X ˆ 1(k - SI|k) as follows:

X ˆ 1(k - SI|k) ) X ˆ 1′(k - SI|k) + K1(k - SI)(y(k - SI) - yˆ (k - SI)) 7. The states X ˆ 1(k - SI|k) are now filtered by EKF1 one sampling instant at a time until the current sampling instant k is reached (X ˆ 1(k|k)). As a result, EKF1 is implemented SI times during this step. Note that EKF1 uses all measurements that are available at each successive sampling instant. The state X ˆ 1(k - SI + 1|k) is stored in a buffer for use at the next sampling instant in step 5. 8. Predict X ˆ 1(k + 1|k) using the process model and calculate P1(k + 1|k), K1(k + 1), and P1(k + 1|k + 1) using the first equation from the sets of equations (4), (5), and (3), respectively. This step also involves the calculation of A(k + 1, k) and H(k + 1). 9. At the next sampling instant, i.e., k ) k + 1, repeat from step 1. In the above algorithm, SI should at least be equal to the largest delay (in this case d) associated with the slow measurements. By doing this, the slow measurements will be used in the reintegration of EKF1 from k - SI to k. If the SI is greater than the measurement delay by an integer n, then the slow measurements will be used n + 1 times by the algorithm. Also, all of the fast measurements are used SI times by the algorithm. The proposed algorithm is illustrated diagrammatically in Figure 1 for SI ) 4 and d ) 2. 3.2. Implementation Issues. 3.2.1. Choice of A and H during Reintegration of Smoothed States. There are three different ways of selecting the matrices A and H in step 7 of the algorithm: (1) use of matrices A(k,k - 1) and H(k) for all sampling intervals between k - SI and k; (2) use of matrices A(k - i, k - i - 1) and H(k - i) ∀ 0 e i e SI, which were calculated at each sampling instant k - i in step 8; (3) the nonlinear model

1040 Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997

is relinearized at each sampling instant k - i ∀ 0 e i e SI, based on the updated value of X ˆ 1, to compute the A and H matrices. The first choice assumes that the process behaves linearly over the smoothing interval period and may not give good results if this assumption is invalid. The second choice is a reasonable compromise between accuracy and computational load. The third choice should result in the best performance, as it relinearizes the process model at each and every sampling interval, based on the most recent state estimates. However, this requires the largest amount of computation. The third option was used in the simulations described in section 5. 3.2.2. Justification and Computation of K2′ i . In step 3, the algorithm adjusts the states of EKF2 by incorporating the slow delayed measurements when they become available. This adjustment directly affects the state vector X ˆ 1(k - SI|k). Any errors in this correction of the states are minimized during reintegration of the state vector up to X ˆ 1(k|k) in step 7 of the algorithm. The computation of K2′ i (k - d) in step 2 is implemented as follows. Whenever a slow measurement is taken, EKF2 calculates the K2i matrix assuming that all measurements were available at that sampling instant. The columns corresponding to the actually available measurements in K2i are used in step 4, and the remaining columns are stored in the buffer (K2′ i ). The model-based output predictions for these delayed measurements (yˆ s) are also stored in a buffer. Whenever these delayed measurements become available, the K2′ ˆ s are used in step 3. i matrix and y 3.2.3. Real-Time Implementation. Steps 1-7 of the proposed algorithm would require computational time on the order of seconds (486DX2/66) between the time the measurements are taken and the calculation of the new state estimates. At this point, any process control algorithm would be executed followed by step 8 during the intersample time. This computational load could be alleviated by modifying the point at which the control algorithm is executed. Once the fast measurements yf(k) are available, step 1 is executed for calculating X ˆ 1(k|k) using the Kalman gain calculated in (the previous) step 8 of the algorithm. At this point, the process control calculations can be performed and the next control move can be implemented immediately based on these estimates. Then, the smoothing and reintegration of the process model (all the remaining steps) are carried out during the intersample time. If the estimator converges to the true state values, the difference between the “ideal” implementation and the suggested modification for real-time applications will be negligible. 3.2.4. State Selection for EKF2. The selection of states for EKF2 is another important issue for this algorithm. Not all states in the original system (eq 1) have to be smoothed. Only a subset of the states associated with EKF1 need to be selected for smoothing, in order to decrease computational load and memory requirements. The set of states in X ˆ 2 must satisfy observability requirements with the measurements being used by EKF2. The main criterion is to include states in EKF2 which are not affected by the fast measurements in order to improve convergence of their estimates. 3.2.5. Multirate Systems. For clarity only, the algorithm was presented for a two-rate measurement

system. This can easily be generalized to multirate systems (measurements taken at different rates and available with different delays). In the proposed algorithm, the delayed measurements are used by EKF1 as they become available during reintegration of the states from sampling interval k - SI to k in step 7. The algorithm does not even require the availability of a delayed measurement at a fixed amount of delay. Therefore, the algorithm lends itself naturally to applications with variable measurement delays. 3.3. Algorithm Simplification. The incorporation of slow measurements in steps 2 and 3 can easily be simplified for certain cases. For strongly observable systems, with large smoothing interval SI, use of a steady-state Kalman filter for EKF2 will not adversely affect the performance of the system. This is due to the fact that EKF2 eventually only affects the states at time k - SI. In fact, if SI is large, K2i will tend to converge to their steady-state values. The smoothed states will then be used by the dynamic EKF1 to reintegrate the process from time k - SI to time k, using all of the available fast and slow measurements. Therefore, for some systems, use of a steady-state smoother EKF2 may not significantly affect the performance of the algorithm but will significantly decrease the complexity and computational load of the algorithm. In the simulations described in section 5, the full algorithm without any simplifications was used. 4. Acrylonitrile/Butadiene Batch Emulsion Copolymerization Process Model The proposed algorithm will now be evaluated on an acrylonitrile/butadiene emulsion copolymerization batch process. The process model characteristics are described below. 4.1. Model. A detailed model for acrylonitrile (AN)/ butadiene (Bd) emulsion copolymerization (NBR or nitrile rubber production) has been developed by Dube et al. (1996). This comprehensive model allows for both thermal and redox initiation (“hot” and “cold” recipes). It accounts for both micellar and homogeneous nucleation, as well as for both water-phase polymerization and reactions inside the polymer particles. In addition, reactions with both water- and monomer-soluble impurities are included. The molecular weight part of the model accounts for transfer to chain-transfer agent, transfer to polymer reactions, and terminal and internal double-bond polymerization in order to explain NBR branching characteristics. The method of moments is used for the calculation of molecular weight averages and average number of branch points. The model consists of 20 nonlinear ordinary differential equations solved simultaneously using a stiff “ode” equation solver. Details about the model development and results from a sensitivity analysis (based on industrial information) can be found in Dube et al. (1996) and Mutha (1996). 4.2. Operating Conditions. The recipe used for the isothermal, AN/Bd emulsion copolymerization batch process is shown in Table 1. The initial values used by the estimator are presented in Table 2. The last entry in Table 1 (monomer-soluble impurities (MSI) flow rate into reactor) represents an impurity disturbance continuously entering the system over the course of the batch which could be caused by an impurity leak, for example. 4.3. Measurements. The simulation results of section 5 have been generated assuming that readily

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1041 Table 1. Recipe Used for Batch Acrylonitrile/Butadiene Emulsion Polymerization ingredient

initial true value

acrylonitrile butadiene water emulsifier initiator activator temp chain-transfer agent MSI no. of particles/L of solvent vol of polymer particles acrylonitrile in particle phase butadiene in particle phase MSI flow rate into reactor

54700.0 mol 114 000.0 mol 906 200.0 mol 675.0 mol 27.0 mol 50.0 mol 285.0 K 100.0 mol 0.0 mol 0.0 0.0 L 0.0 mol 0.0 mol 6.0 × 10-5 mol/min

Table 2. Covariance Matrices for Fixed-Lag Smoothing-Based EKF state

initial variance (P)

state variance (Q)

state initialization

NCTA NMSI PA PB 10-16Np Vp Q0 Xs1 Xs2

Covariance State Matrices for EKF1 1.0 × 10-1 1.0 × 10-2 100 mol 1.0 × 10-7 1.0 × 10-9 0 mol 1.0 × 101 1.0 × 10-2 0 mol 2 -1 5.0 × 10 1.0 × 10 0 mol 1.0 × 10-2 1.0 × 10-7 0 1.0 × 10-1 1.0 × 10-2 0L 1.0 × 102 1.0 × 10-1 0 mol 1.0 × 10-5 1.0 × 10-8 1 5.0 × 10-7 1.0 × 10-8 1

PB 10-16Np Q0

Covariance State Matrices for EKF2 1.0 × 102 1.0 × 10-1 1.0 × 10-2 5.0 × 10-4 1.0 × 102 5.0 × 10-4 measurement

variance (R)

density (F) particle diameter (Dp) PA ln viscosity

1.0 × 10-7 1.0 × 10-10 1.0 × 10-4 5.0 × 10-5

available sensors (see Chien and Penlidis (1990)) are used. The density (F) of the polymer and, hence, conversion can be measured instantaneously and are therefore used as a fast measurement in the simulations. Percentage-bound acrylonitrile (PA) in the polymer (a measure of copolymer composition), average particle size of the particles (Dp), and a measure of the polymer viscosity (ν) are all slow delayed measurements. The polymer viscosity is related to the number-average molecular weight by a nonlinear correlation. In the simulations, the base sampling interval was chosen as 5 min, with density available every sampling interval without delay. The sampling interval for the slow measurements and the measurement delay associated with these measurements ranged from 6 to 30 times the base sampling interval. 4.4. Estimator States and Observability. The process states selected for estimation were based on observability conditions. The states selected for the state estimator EKF1 are the moles of CTA, MSI, [MA]p (or PA), [MB]p (or PB), Np, Vp, and Q0. Two nonstationary states are also considered for inclusion in the estimator to achieve integral action in the state estimates (Stephanopoulos and San, 1984; Kozub and MacGregor, 1992). One of these nonstationary states (Xs1) was multiplicative with the total number of particles (Np), and the second of these states (Xs2) was multiplicative with the volume of polymer particles (Vp). The states affecting conversion are observable by using the fast measurement (density). However, esti-

mator states like Np and Q0 (zeroth moment of the chain length distribution) do not have any effect on the fast measurement (density) and, therefore, will not be observable using the fast measurement alone. During the filtering part of the algorithm when only fast measurements are used, these unobservable states will naturally have a very low Kalman gain and, therefore, are almost unaffected. Since the state Q0, which directly affects the number-average molecular weight, and Np are not observable from the density, these states will be included in EKF2. In addition, percentagebound butadiene (PB) has also been included in EKF2. In this polymerization system, the weight-average molecular weight and trifunctional and tetrafunctional branching frequencies are not observable. Including an appropriate sensor (in the future) could make some of these states observable. However, in these simulations, any states not included in the estimator will be predicted solely using the model equations. Measurement signals are often noisy. Also, unaccountable deviations from ideal conditions occur in almost every application. These effects can be simulated by the addition of noise in the measurements and process states. Independent white noise sequences were added to all the measurements and to the state equations for I, S, MSI, PA, and PB. These states were chosen as they affect all of the important product properties. The noise sequences were designed by multiplying a white noise sequence of unit variance by a percentage (σ) of the corresponding state or measurement. 5. Simulation Results The simulation examples illustrate the features of the proposed algorithm. They also show the algorithm’s ability to effectively work with the slow measurements currently available in the polymer industry (Chien and Penlidis, 1990). The nonstationary state Xs1 was always included, but unless otherwise specified, the nonstationary state Xs2 was not included. The simulation plots that follow present the most important variables of the system. MSI in the system is not observable with the considered measurement vector and, therefore, never converges to its true value. However, the main effect of MSI is on the number of particles generated, which is observable. The estimator accounts for the effect of MSI by manipulating the number of particles and other related parameters of the system. The tuning matrices and the initialization of estimator states are listed in Table 2. Tuning of the estimator was carried out in a time-consuming, trial-and-error manner. 5.1. Base Case. Simulation runs were performed with all three slow measurements being taken every sixth sampling interval. The measurement delay associated with these measurements was also set as six sampling intervals (30 min). State and measurement noise at 6% ()3σ) of the respective states and measurements was used. A smoothing period of 12 sampling intervals was selected. The base case simulation was performed with a flow rate of 6.0 × 10-5 mol/min of MSI into the reactor. The performance of this base case is contrasted with the case where there is no monomersoluble impurity being added to the process, i.e., an MSI flow rate of 0 mol/min. The presence of MSI has an often ignored and unappreciated effect on particle nucleation and, consequently, on other polymerization variables (Huo et al., 1988; Penlidis et al., 1988). The performance of these runs is contrasted in Figure 2.

1042 Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997

Figure 2. Base case and effect of impurity on convergence of the estimator. MSI ) 6.0 × 10-5 (estimator (- ‚ -), true (s)) (base case) and MSI ) 0.0 mol/min (estimator (•), true (- -)).

The number of particles generated in these cases is substantially different. In both cases, the estimator states converge to their true values, except for the particle diameter. The convergence of the estimator for the particle diameter will be discussed in more detail in section 5.4. The results in Figure 2 illustrate the fact that the estimator effectively handles the presence of MSI, and in general, the results agree with experimental results found in the literature (Huo et al., 1988; Penlidis et al., 1988). 5.2. Effect of Smoothing Interval. Here, the base case conditions were used except with state and measurement noise at 15% ()3σ) of the respective states and measurements. The smoothing interval was varied in the simulations. Figure 3 shows the estimator performance with smoothing intervals of 6, 7, and 12. The effect of using a larger SI is most visible on the estimated number-average molecular weight where convergence to the true state is faster with smaller oscillations at a larger SI. Beyond a particular level of smoothing (7 in this case), the performance improvement is negligible. 5.3. Effect of Measurement and State Noise. The estimator performance was evaluated in this case at an

Figure 3. Effect of smoothing interval on convergence of the estimator. Smoothing interval of 6 (O), 7 (+), and 12 (×) against true (s) values.

even higher noise level than used previously. The run presented in section 5.2 was with 15% (3σ) noise in both

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1043

Figure 4. Robustness to noise level. Measurement and state noise of 30% (3σ). Smoothing interval of 6 (O) against true (s) values.

Figure 5. Robustness to noise level. Measurement and state noise of 30% (3σ). Smoothing intervals of 12 (O) and 18 (+) against true (s) values.

measurements and estimator states. These runs have been repeated here with 30% (3σ) noise in both measurements and estimator states. The estimator diverges and exhibits numerical difficulties for a smoothing interval of 6, after 50 sampling intervals (Figure 4). With smoothing intervals of 12 and 18 (Figure 5), it is clear that the estimator’s performance improves with increasing smoothing interval. At this noise level, it was found that the performance of the estimator was marginal with a smoothing interval of 7. This set of simulations shows that the use of larger smoothing intervals results in better rejection of state and measurement noise by the estimator. 5.4. Effect of Nonstationary States. In all the simulation results presented in sections 5.2 and 5.3, an offset was observable between the true and estimated values of the particle diameter. In this case, it is demonstrated that this offset can be reduced through the addition of Xs2. The effect of nonstationary states has been discussed by Stephanopoulos and San (1984) and Kozub and MacGregor (1992) and experimentally verified by Dimitratos et al. (1991). Here, the base case conditions were used with measurement and state noise at 30% ()3σ) of the respective states and measurements. The performance of this case is compared with the case

Figure 6. Effect of including an additional nonstationary state. Smoothing interval of 12, with (+) and without Xs2 (O) against true (s) values.

where Xs2 was not included in Figure 6. The average particle diameter shows better convergence to its true value with the additional nonstationary state. 5.5. Effect of Model-Process Mismatch. In practice, model-process structural mismatch and errors in the values of some of the parameters of the model are to be expected. These mismatches should be compensated for by the estimation algorithm. In this case, measurement and state noise of 15% ()3σ) was used, and both nonstationary states were included. Simulation runs were performed with all three slow measurements being taken every 16th sampling interval. The measurement delay associated with these measurements was also set as 16 sampling intervals (80 min). Smoothing intervals of 16, 24, and 32 sampling intervals were examined. All runs were performed with a MSI level corresponding to 6.0 × 10-5 mol/min (same as base case). Changes in the propagation rate constants are known to significantly affect the reaction kinetics. In the true process model, the lumped propagation rate constant was increased by a factor of 30% compared to the base case. The performance of these runs is shown in Figure 7. The estimator performance with respect to the numberaverage molecular weight is a function of the smoothing interval. The estimator does not converge for a smoothing interval of 16 and has rather poor convergence with a smoothing interval of 24. However, it does converge with a smoothing interval of 32. As in all previous simulations, good convergence is observed for smoothing intervals greater than the largest sum of the sampling interval plus measurement delay for each slow measurement considered separately. This is because when this choice for SI is made, at least one set of slow measurements is used at each and every sampling interval by the estimator. For the simulations presented in this case, the largest sum is equal to 16 + 16 ) 32. Whenever there is model-process mismatch, the estimator will try to minimize the effect of the mismatch by biasing some of the states. The states which become (artificially) biased depend on the set of measurements and the set of states being used in the filter. A large part of this bias is normally absorbed by the nonstationary states. In these simulations, the number of particles in the reactor has been multiplied by a

1044 Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997

Figure 7. Effect of model-process mismatch in propagation rate constant by 30%. Smoothing interval of 16 (O), 24 (+), and 32 (×) against true (s) values.

nonstationary state (Xs1) and has, therefore, converged to a biased value. The performance of the estimator for a system with model-process mismatch can be contrasted with the base case. For both the base case (Figure 2) and the simulation with model-process mismatch (Figure 7), the estimator effectively tracks the number-average molecular weight (M h n). However, this is achieved in very different ways. M h n is the ratio of Q1 to Q0, multiplied with the effective molecular weight, Mweff. The Mweff is calculated as the sum of the individual monomer molecular weights weighted by the copolymer composition and can be directly related to the mole fraction reacted acrylonitrile (FA) as follows:

Mweff ) MWmAFA + MWmB(1 - FA) The Q0 and Q1 for both cases are shown in Figure 8. In the base case, both Q0 and Q1 converge to their true values. However, when there is a mismatch in the propagation rate constant, they do not converge to their true values. Q0 is a state in the estimator, whereas Q1 is not. When the propagation rate constant is increased, M h n increases. Since Q1 is not a state in the estimator, the estimate of Q1 based solely on the model equations remains unchanged, even though the kinetics of the true process have changed. The estimator responds by biasing Q0 and FA such that convergence of M h n is achieved. 5.6. Effect of Very Slow Measurements. This case evaluates the effect of having very slow measurements on the convergence of the estimator states. Runs were performed with slow measurements being taken after 30 sampling intervals. The measurement delay was also set equal to 30 sampling intervals. State and measurement noise at 15% ()3σ) of the respective states and measurements was added, and both nonstationary states were included. Smoothing intervals of 30, 45, and

Figure 8. Performance of the estimator in tracking Q0 (estimated (O), true (s)) and Q1 (estimated (+), true (- -)). (a) Base case and (b) mismatch in propagation rate constant.

60 were used by the estimator. The same 30% mismatch in the lumped propagation rate constant was present in these simulations. The performance of the estimator for smoothing intervals of 30, 45, and 60 is shown in Figure 9. The estimator is not able to achieve good convergence with smoothing intervals of 30 and 45 for the numberaverage molecular weight. The estimator diverges for the smoothing interval of 60 in the number-average molecular weight. It is worth noting here that without the propagation rate constant mismatch, the estimator had much better convergence properties for all SI values. By taking more frequent measurements but without requiring any decrease in the measurement delay itself, convergence can be greatly improved even in this difficult case. Slow measurements taken every 30 sampling intervals (every 150 min) having a measurement delay of 30 sampling intervals (150 min) imply that only one set of measuring devices is being used to analyze all of the samples in series. If two sets of measuring devices were available, the sampling frequency could be doubled. For this case, it would mean that a measurement could be taken every 15 sampling

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1045 Table 3. Covariance Matrices for the Basic EKFa state

initial variance (P)

state variance (Q)

NCTA NMSI PA PB 10-16Np Vp Q0 Xs1

1.0 × 10-1 1.0 × 10-7 1.0 × 101 5.0 × 102 1.0 × 10-2 1.0 1.0 × 102 1.0 × 10-4

1.0 × 10-4 1.0 × 10-9 1.0 × 101 5.0 × 10-1 1.0 × 10-7 1.0 × 10-2 1.0 1.0 × 10-11

a Same measurement variances as fixed-lag smoothing-based EKF.

Figure 9. Effect of very slow measurements. True values (s). (a) Smoothing intervals of 30 (O) and 45 (+). (b) Smoothing interval of 60 (O). (c) Using two sets of measuring devices with a smoothing interval of 45 (O).

intervals (75 min), although the measurement will not be available for 30 sampling intervals (150 min). The performance of the estimator with slow measurements taken every 15 sampling intervals and a measurement delay of 30 sampling intervals using a smoothing interval of 45 is also shown in Figure 9. Now, the estimate of M h n converges. These simulations demonstrate that using multiple sets of measuring devices for measurements with large measurement delays can greatly improve on-line state estimation. This situation is often seen in industry where some key measurements only become available after a long period of laboratory analysis and are then considered obsolete (Ogunnaike and Gopalratnam, 1991). However, the above results suggest an alternative where, in spite of the long analysis time, the measurements can effectively be used for achieving good results with the proposed estimator. 5.7. Comparison with Standard EKF. The performance of the proposed estimator was finally compared with the standard EKF. The standard EKF here refers to the filter implementation without any smoothing. Since this algorithm cannot cope with the delayed measurements, it will operate with the multirate measurements assuming all measurements are available without measurement delay. This means that in the comparisons to follow, the standard EKF has a significant advantage over the proposed algorithm in terms of availability of measurements. The simulations with the standard EKF were performed with the same conditions as the base case described in section 5.1. The standard EKF was tuned separately in order to achieve its “best” performance. The tuning matrices are listed in Table 3. The performance of the standard EKF has been compared with that of the proposed algorithm for the base case conditions in Figure 10. Although the standard EKF has the

advantage of working with undelayed measurements, its performance is poorer. Much better performance is achieved using the fixed-lag smoothing-based EKF since one of its main advantages is the ability to reuse a measurement multiple times. The algorithm by Ellis et al. (1988) would at best approach the performance of the standard EKF with all the measurements available without measurement delay. Therefore, the simulation results of Figure 10 illustrate the advantage of using the proposed algorithm over that of Ellis et al. (1988). 5.8. Features of the Proposed Algorithm. Several important features of the proposed fixed-lag reiterated EKF have been illustrated through an extensive case study, and the results have important implications for industrial practice. First, the estimator exhibits generally good convergence to the true process states and robustness to measurement and process noise. Second, the estimator is capable of handling systems with multirate measurements. For systems where measurement delays are very long and may be expected to cause estimator divergence, the simulations demonstrate that estimator convergence can be achieved by taking more frequent measurements and using multiple sets of measurement devices for their analysis. Third, the proposed algorithm can accept measurements with a variable measurement delay. This is because slow measurements are used only during the reintegration of the smoothed states. At a given sampling instant, if a measurement is available, the algorithm will use it for filtering. If it is not available, the filtering will still be performed with the other measurements available. Of course, there is a cost associated with obtaining these improvements in that the proposed algorithm requires a rather large computational effort and this computational load increases linearly with the smoothing interval. The memory requirements for the algorithm are higher than the standard EKF as all of the measurements, selected states, and associated Kalman gain matrices are stored over a moving time window. However, these costs are not significant even at this time with currently available computers. The other design cost with using the proposed algorithm is the need to select the states for the smoothing part of the algorithm and to tune the smoother. Through usually not mentioned, the implementation of such complex algorithms through programming in high-level languages takes considerable time and effort. Selection and use of proper tools, techniques, and programming languages are all important factors. The proposed algorithm was predominantly written in C in order to use its flexible dynamic allocation features. The programs were designed to be very flexible in terms of changing smoothing intervals, adding/deleting process states, changing the interval at which slow measurements are taken and their associated measurement

1046 Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997

Figure 10. Comparison of the performance of basic EKF (- -) algorithm with the base case (•) of the proposed algorithm against true (s) values.

delays, etc. Use of an object-oriented programming language like C++ would be a good choice for implementing the proposed algorithm. One final note should be made about the fact that the algorithm performance has been illustrated through simulations on a specific process. Nonetheless, the emulsion copolymerization multirate process represents a very challenging and realistic example, and similar performance can be expected with other systems. 6. Concluding Remarks A fixed-lag smoothing-based extended Kalman filter algorithm was proposed for systems with multirate measurements. The performance of the filter has been evaluated through a case study on a complex emulsion copolymerization batch process. The proposed estimator uses measurements multiple times to achieve better convergence and is robust to state and measurement noise. The estimator can also handle systems with large measurement delays. For cases where the estimator diverges due to large measurement delays, the simulations demonstrate that convergence can be achieved by using multiple measurement devices. The simulations

also demonstrate that the performance of the proposed algorithm using delayed measurements is superior to the standard extended Kalman filter with all measurements being available instantaneously. Acknowledgment The authors acknowledge financial support from the Natural Sciences and Engineering Research Council of Canada and also thank Uniroyal Chemical, USA, and ICI/Glidden, Worldwide, for useful discussions. Nomenclature (‚ˆ) ) estimate of (‚) (‚)1 ) (‚) parameter related with EKF1 (‚)2 ) (‚) parameter related with EKF2 A ) acrylonitrile B ) butadiene CTA ) chain-transfer agent Dp ) diameter of particles (dm) FA ) mole fraction reacted acrylonitrile Ki ) Kalman gain matrix for ith state(s) m ) monomer phase (as a subscript) M h n ) number-average molecular weight

Ind. Eng. Chem. Res., Vol. 36, No. 4, 1997 1047 MSI ) monomer-soluble impurities Mweff ) effective molecular weight Ni ) number of moles of species i n j ) average number of radicals per particle Np ) total number of particles in the reactor (per liter of water) Pi ) moles of species i in particle Pij ) covariance of ith states(s) with jth state(s) Q0 ) zeroth moment of chain length distribution (mol) Q1 ) first moment of chain length distribution (mol) Vp ) total volume of polymer particles (L) X ˆ (k) ) state estimates at time k Xs ) nonstationary state X ˆ (k|p) ) estimate of state x at sampling interval k based on measurements until sampling interval p yf ) fast (frequent) measurements ys ) slow (infrequent) measurements F ) density (g/L) ν ) polymer viscosity (P) σ ) standard deviation

Literature Cited Adebekun, D. K.; Schork, F. J. Continuous solution polymerization reactor control. 2. Estimation and nonlinear reference control during methyl methacrylate polymerization. Ind. Eng. Chem. Res. 1989, 28, 1846-1861. Anderson, B. D. O.; Moore, J. B. Optimal Filtering; Prentice-Hall: Englewood Cliffs, NJ, 1979. Chien, D. C. H.; Penlidis, A. On-line sensors for polymerization reactors. J. Macromol. Chem. Phys. Rev. 1990, C30, 1-42. Dimitratos, J.; Georgakis, C.; El-Aasser, M. S.; Klein, A. Dynamic modeling and state estimation for an emulsion copolymerization reactor. Comput. Chem. Eng. 1989, 13, 21-33. 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, 3203-3218. Dimitratos, J.; Elicabe, G.; Georgakis, C. Control of emulsion polymerization reactors. AIChE J. 1994, 40, 1993-2021. Dube, M. A.; Penlidis, A.; Mutha, R. K.; Cluett, W. R. Mathematical modeling of emulsion copolymerization of acrylonitrile/ butadiene. Ind. Eng. Chem. Res. 1996, 35, 4434-4448. Elbert, T. F. Estimation and Control of Systems; Van Nostrand Reinhold: New York, 1984.

Elicabe, G. E.; Meira, G. R. Estimation and control in polymerization reactors. Polym. Eng. Sci. 1988, 28, 121-135. Ellis, M. F.; Taylor, T. W.; Gonzalez, V.; Jensen, K. F. Estimation of the molecular weight distribution in batch polymerization. AIChE J. 1988, 34, 1341-1353. Gudi, R. D.; Shah, S. L.; Gray, M. R. Adaptive multirate state and parameter estimation strategies with application to a bioreactor. AIChE J. 1995, 41, 2451-2464. Huo, B. P.; Campbell, J. D.; Penlidis, A.; MacGregor, J. F.; Hamielec, A. E. Effect of impurities on emulsion polymerization: case II kinetics. J. Appl. Polym. Sci. 1988, 35, 20092021. Kozub, D. J.; MacGregor, J. F. State estimation for semi-batch polymerization reactors. Chem. Eng. Sci. 1992, 47, 1047-1062. MacGregor, J. F.; Penlidis, A.; Hamielec, A. E. Control of polymerization reactors: a review. Polym. Process. Eng. 1984, 2, 179-206. Meditch, J. S. A survey of data smoothing for linear and nonlinear dynamic systems. Automatica 1973, 9, 151-162. Mutha, R. K. Nonlinear estimation and control: Applications to polymer reactors. Ph.D. Thesis, University of Toronto, 1996. Ogunnaike, B. A.; Gopalratnam, P. C. A two-tier multirate control system structure for an industrial terpolymerization reactor. Annual AIChE Conference, Los Angeles, 1991, Paper 145h. Penlidis, A.; MacGregor, J. F.; Hamielec, A. E. Effect of impurities on emulsion polymerization: case I kinetics. J. Appl. Polym. Sci. 1988, 35, 2023-2038. Schuler, H.; Suzhen, S. Real-time estimation of the chain length distribution in a polymerization reactor. Chem. Eng. Sci. 1985, 40, 1891-1904. Schuler, H.; Schmidt, C. U. Model-based measurement techniques in chemical reactor applications. Int. Chem. Eng. 1993, 33, 559-576. Stephanopoulos, G.; San, K. Y. Studies on on-line bioreactor identification. Biotech. Bioeng. 1984, 26, 1176-1188.

Received for review February 20, 1996 Revised manuscript received July 2, 1996 Accepted July 2, 1996X IE9601007

X Abstract published in Advance ACS Abstracts, February 15, 1997.