Ind. Eng. Chem. Res. 2004, 43, 7227-7237
7227
Regulation of the Emulsion Particle Size Distribution to an Optimal Trajectory Using Partial Least Squares Model-Based Predictive Control Myung-June Park, Mustafa T. Dokucu, and Francis J. Doyle III* Department of Chemical Engineering, University of California, Santa Barbara, Santa Barbara, California 93106
This paper describes a statistical model-based control strategy for the regulation of the particle size distribution in emulsion copolymerization. A database is generated using both a mechanistic model and experimental data. A partial least squares method is applied at each sample time to generate the prediction model to be used in a predictive control algorithm. The proposed strategy is implemented in an experimental reactor in which the copolymerization of vinyl acetate/butyl acrylate is carried out. The controller is effective in updating the input trajectories to regulate the controlled distribution to a desired target in the presence of strong disturbances. 1. Introduction The population balance equation (PBE) has been widely used in the development of mechanistic models for chemical processes including emulsion polymerization, aerosols, granulation, crystallization, and precipitation (see work by Ramkrishna1 for the theoretical aspects of the PBE). The resulting models have been used for process optimization and, recently, for modelbased control of the underlying distribution (e.g., size, shape, etc.). One limitation in the application of PBEs to the control of distributions is that the discretized (lumped) versions of the model are typically high-order, ill-conditioned, and uncontrollable. In the discretized version of the PBE, the differences of the order of magnitude between elements are so severe that the Jacobians become hypersensitive to slight changes in the variables and, thus, optimization based on the gradient method is not applicable for the calculation of the optimal input values. Furthermore, because the magnitude changes very rapidly, it is difficult to find the bounds to normalize the variables.2 Consequently, most of the published work on the control of the distributed parameter systems have focused on moment approximations. Kalani and Christofides3 applied a moment method in the prediction of the particle size in an aerosol process and controlled the geometric average particle diameter using a Luenberger-type observer combined with an output feedback controller. Zeaiter et al.4 used the particle size polydispersity index in an objective function in order to calculate the optimal trajectory of the monomer feed rate. Valappil and Georgakis5 introduced a correlation between the properties of the particle size and average molecular weight and the end-use properties of the tensile strength and melt index, and a nonlinear model predictive control based on successive target region linearization was applied to the control of these end-use properties. To overcome the computational limitations of the mechanistic model, chemometrics-based models have been suggested as an alternative.6 Among the methods * To whom correspondence should be addressed. Tel.: +1 (805) 893-8133. Fax: +1 (805) 893-4731. E-mail: doyle@ engineering.ucsb.edu.
based on chemometrics, the partial least squares (PLS)7 model has been most popularly used as an alternative to the more classical multiple-linear regression (MLR) and principal-component regression (PCR) methods, because this method deals with the problems of collinearity and the high dimension of data blocks. This characteristic is particularly appropriate for chemical processes because chemical processes are intrinsically multivariable and consist of mutually correlated variables. Standard linear PLS involves the projection of X [which relates to the process variables (inputs)] and Y [which corresponds to the quality or reference variables (outputs)] onto a subset of latent variables that are referred to as the input and output scores. Dayal and MacGregor8 modified the original kernel algorithm9 to generate a more efficient algorithm, while other researchers10,11 used a singular value decomposition to compute components of the PLS kernel without using the time-consuming iterative procedure of nonlinear iterative PLS. As discussed by Baffi et al.,12 when applying linear PLS to a nonlinear problem, the minor latent variables cannot necessarily be discarded because they may describe not only the noise in the data but also significant components of the variance-covariance structure. Identifying the relevant latent variables to incorporate into the model requires additional analysis, and the final model may contain too many components. Therefore, a nonlinear variant of the PLS has been developed. A simple approach is to extend the input matrix by including nonlinear combinations of the original variables such as logarithms,10 square roots, and cross products.13 Durand11 and Baffi et al.12 carried out a nonlinear transformation of the input data using B splines.14 The problem with this method is that when the number of variables increases, the number of square roots or cross products rapidly increases. Recently, Wold et al.15 described an orthogonal signal correction (OSC) method in which scores orthogonal to Y are found to eliminate Y-irrelevant data from X, and they showed that it is always possible to find an exactly orthogonal OSC solution when the number of X variables is larger than the number of training samples. Trygg and Wold16 proposed an improved preprocessing method, called
10.1021/ie034304g CCC: $27.50 © 2004 American Chemical Society Published on Web 05/15/2004
7228
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
orthogonal projections to latent structures (O-PLS) obviating an iterative search for scores orthogonal to Y. Another approach for the nonlinear PLS approach involves a nonlinear relationship in an inner mapping between input and output matrices. Wold et al.17 used a polynomial while Holcomb and Morari18 applied a neural network for the nonlinear inner relationship. PLS model-based control has been widely used in chemical process applications. Elizalde et al.19 applied PLS methods to construct a model relating the adhesive properties to the full molecular weight distribution (MWD), and they demonstrated the production of polymers with the desired MWD and copolymer composition using the feed profiles of monomers and chain-transfer agents as manipulated variables. Flores-Cerillo and MacGregor13 controlled the full particle size distribution (PSD) of a styrene emulsion polymerization reactor using a batchwise-adaptive PLS model by which injections of emulsifier at different times were selected, while Doyle et al.20 combined a mechanistic PBE model with a PLS model to compensate for uncertainty in the coagulation mechanism as well as more general model mismatch. This hybrid model was utilized to calculate the optimal trajectories of surfactant and initiator using a successive quadratic program for the control of the PSD. In a batch pulp digester, Kesavan et al.21 compared the performance of a controller in which the optimal inputs of the cooking temperature and batch end time were calculated directly using the PLS model with that of a nonlinear H-factor model-based controller. In this study, a PLS method is utilized to predict the weight-averaged PSD at the final time of a batch as a function of the future input trajectory. The prediction models based on the PLS method are obtained at each sample time and are used to calculate optimal trajectories using a predictive control algorithm. The PLS model-based predictive control is introduced in section 2, and the mechanistic model of the semibatch vinyl acetate (VAc)/butyl acrylate (BuA) emulsion copolymerization system is summarized in section 3. The details of the experimental system and online measurement are given in section 4, followed by the results and discussion in section 5. 2. PLS Model-Based Predictive Control A general method for construction of a model with the relationship Y ) f(X) is multiple-linear regression (MLR).7 This can be represented mathematically as
Y ) XB + EMLR
(1)
where X, Y, B, and EMLR represent the regressor (or input) block, the response (or output) block, the regression (or sensitivity) matrix, and the residual matrix, respectively. Depending on the number of samples (n) and the number of independent variables (m), one can distinguish three possible cases. First, when m > n, there are an infinite number of solutions for the regression matrix (B), while in the case of m ) n, a unique solution for B is provided (though this is rare in practice). In the third case, where m < n, there is no exact solution, but an approximate solution can be determined by minimizing the magnitude of the residual vector (EMLR). The leastsquares solution is given by
B ) (XTX)-1XTY
(2)
This equation reveals a common problem in MLR; namely, the inverse of XTX may be difficult to compute if the input data block is of high dimension and its elements are highly correlated, which is denoted collinearity or singularity. The PLS approach provides an alternative to overcome this problem. PLS methods involve an approximation of the X and Y spaces with their respective score matrices and the maximization of the correlation between the original data blocks.6,7 A simplified model consists of a regression between the scores for the X and Y blocks. The PLS model can be interpreted as consisting of outer relations (X and Y blocks individually) and an inner relation (linking both blocks). The outer relation for the X block is a
T
X ) TP + E )
thphT + E ∑ h)1
(3)
One can build the outer relation for the Y block in the same manner: a
T
Y ) UQ + F* )
uhqhT + F* ∑ h)1
(4)
Here, T and U represent the score matrices that explain the coordinates of the respective points on the principal components, while P and Q denote the loading matrices that determine the directions of the principal components. E and F* are residual matrices. Each pair of latent variables accounts for a certain amount of the variability in the input and output blocks. The first a latent variables account for most of the variance of the data blocks, and the other latent variables capture measurement and/or process noise in the data. It is worth noting that the latent variables corresponding to small eigenvalues should be left out in order to avoid collinearity problems.22 In the limit of enough latent variables, E ) F* ) 0. However, it is the intention to both describe Y accurately (make |F*| as small) and, at the same time, construct an efficient representation between X and Y. A simple criterion for model building is to use a threshold value for F*, but cross validation can be used to determine the number of latent variables for a robust model.12 Cross validation6 is performed by an iterative procedure in which a fixed fraction of the elements are first removed from the data matrix; typically in some cases, it is preferable to delete diagonal lines of elements. When the algorithm is initialized with one principal component, a principal-component analysis is performed on the reduced data set to extract a components. Then, the deleted data values are predicted using these principal components, and the predictive residual error sum of squares (PRESS) is computed and stored for later use:
PRESS(a) )
∑i ∑j eij2
(5)
where i and j denote the indices for variable and training sets, respectively. As a (the number of components) approaches the true number of significant components, the prediction should improve and, thus, the PRESS should decrease. However, as the truly significant number is exceeded, noise is included with the model, and the predictive ability is reduced. As a result,
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7229
the PRESS increases, and an optimal value of a can be determined. The simplest model for the inner relation is a linear one:
u ˆ h ) bhth
(6)
ˆ hTth/thTth. bh plays the role of the regression where bh ) u coefficient in the MLR and PCR models. The input matrix (X) includes online measurements that are available at the current time (xcurrentT), control inputs at previous sample times (upastT), and inputs to be made at current and future sample times (ufutureT). Flores-Cerillo and MacGregor13 suggested the use of measurements carried out offline (xoffT) and any relevant information from immediately previous batches (xpriorT) for better prediction of the behavior of the current batch. However, xoffT and xpriorT are not used in this study. The output matrix (Y) contains the observations of quality variables (yT) to be controlled. Finally, the prediction model for quality variables for the current batch is constructed by applying the PLS algorithm at each sample time in linear regression form as follows:
ˆi yˆ T ) xiTb
(7)
Here, subscript i denotes the current sample time and circumflex indicates the estimated value. The current state is
xiT ) [xcurrentT, upastT, ufutureT]
(8)
The prediction model is utilized in a control algorithm to calculate optimal future input trajectories. Because the prediction model is a linear function of inputs, the application is straightforward. One can apply the minimum variance control algorithm or linear quadratic regulator.13 In this study, a model predictive control algorithm is used:
ufuture ) arg min{[ysp - yˆ (ufuture)]TΛy × [ysp -yˆ (ufuture)] + ufutureTΛuufuture} (9) subject to C‚ufuture e c
(10)
The proposed control strategy involves the recalculation of the prediction model as new online measurements are available. Hence, the algorithm is based on multiple local models that cover operating regions around the current state. Flores-Cerillo23 suggested the utilization of a global PLS model that covers the entire operating window by constructing the input block including all of the measurements available during the entire batch time. When this model is applied to the computation of future control inputs, a missing data algorithm24,25 is used to estimate the unknown future measurements. However, given the transient and highly nonlinear behavior of batch and semibatch processes, the use of a single global model may lead to a complex prediction model with deteriorated prediction ability. Therefore, a multiple local PLS model-based algorithm13 is used in this study for the control of properties in a semibatch emulsion copolymerization reactor. Johansen and Foss26 suggested an operating regime based process modeling scheme, but deterioration of the performance is expected when the current operating region is close
to the boundaries of prespecified operating regimes. Meanwhile, the local models in this study are dependent on the current state. The proposed algorithm is similar to the model-on-demand predictive controller (MoDPC)27,28 in which the data within a domain of a prespecified distance from the current state are used to build a local model. However, the MoD-PC suffers from the determination of an adequate distance function and the local model structure for a good performance, and it requires a great number of training batches for a successful local model. Furthermore, the MoD-PC cannot be applied to the system in which the measurement is infrequently available. 3. Modeling of a VAc/BuA Emulsion Copolymerization Reactor For the generation of the database for the PLS model, input and output data that cover the operating range are required. However, considering the application of the method in an industrial context in which the operation is conducted within a narrow range, it is difficult to obtain fully excited data. Therefore, the input and output data blocks are generated by a first principles model. Immanuel et al.29 developed a mechanistic model that describes the dynamic behavior and the properties of a VAc/BuA emulsion copolymerization reactor and presented a computationally efficient method to solve the distributed parameter system.30 That model is utilized in this work, and the salient features are outlined here. The particle density function [F(r,t)] is calculated using a population balance model composed of three underlying processes in an emulsion polymerization reactor: nucleation, growth, and coagulation.
∂ ∂ F(r,t) + [F(r,t) Rgrowth(r,t)] ) ∂t ∂r Rnuc(r,t) + Rcoag(r,t) (11) The second term on the left-hand side, the partial derivative with respect to r, is the growth kernel, whose expression is size-dependent in this model. As for the nucleation rate [Rnuc(r,t)], two types of mechanisms are considered. One is micellar nucleation, which occurs when the concentration of free surfactant (Sw) in the aqueous phase exceeds the critical micelle concentration, and the other is homogeneous nucleation, which forces the oligomers in the aqueous phase with length over the critical chain length (jcr) to aggregate surfactant molecules around them. The coagulation rate [Rcoag(r,t)] is calculated from considerations of the thermodynamic instability of colloidal particles. The mathematical forms of the kernels that account for these processes are as follows:
Growth Rgrowth(r,t) )
dr dt
)
1
2
2
∑ ∑kpijpi i)1 j)1
4πr2Fp
n j (r,t) NA
[Mj]pMWj (12)
Nucleation Rnuc(r,t) ) Rmicellar + Rhomogeneous
(13)
7230
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
Figure 1. Schematic diagram of experimental system for VAc/BuA emulsion copolymerization. jcr-1 2
Rmicellar )
∑ ∑ l)0 i)1
l ei,micelle pwi[Pw]lCmicelleVaq
jcr-1 ]Vaq Rhomogeneous ) kw pav[Pw
(14) (15)
Coagulation Rcoag(r,t) )
250 3 sity function as wPSD(ri,t) ) ri3F(ri,t)/∑i)1 ri F(ri,t), it also is composed of 250 elements. Mass balance equations for the initiator, initiator radical, and monomers are considered. On the basis of a redox initiation system employed in the VAc/BuA copolymer system, the following material balances for the oxidizer (Iw), the reducer (Y2), and the initiator radical (Rw) are calculated:
H(rupper-r) Rformation(r,t) - H(rcutoff-r) Rdepletion(r,t)
d([Iw]Vaq) ) -kd1[Iw][Yr1] + υIw dt
(20)
d([Y2]Vaq) ) -rIkd1[Iw][Yr1] + υY2 dt
(21)
(16) Rformation(r,t) ) 1 Vaq
∫rr
r2 β(r′,r′′) F(r′,t) F(r′′,t) 3 dr′ (17) (r - r′3)2/3
max
nuc
Rdepletion(r,t) )
1 Vaq
∫r
d([Rw]Vaq) dt
rmax nuc
β(r,r′) F(r,t) F(r′,t) dr′
(18)
4πD0(r + r′) WS
(19)
β(r,r′) ) c1
The PBE can be solved by a variety of numerical methods including weighted residuals or finite differences. In the present study, the numerical solution method developed by Immanuel and Doyle30 is applied. The continuous PSD is divided into a finite number of sections within which an integral quantity of the distribution is defined, in such a way that nucleation, growth, and coagulation processes that occur in each of the sections are considered individually to update the particle count. Because this method determines the nucleation, growth, and coagulation processes by the underlying thermodynamic and kinetic events, governed by relatively simple equations with relaxed stiffness characteristics, the associated stiffness in the full solution is removed. To discretize the PBE for the particle density function (cf. eq 11), 250 elements (or grids) with a width of 2 nm are used to cover the range of particle sizes spanning from 2 to 500 nm and mass balance equations for the two monomers, aqueous-phase volume, surfactant, oxidizer, reducer, and initial radicals are calculated. Because the weight-averaged PSD (the controlled output) is calculated using the particle den-
2
) kd1[Iw][Yr1] - Vaq
kri[Rw][Mi]w ∑ i)1 jcr-1
w [Rw]( Vaqktav
[Pw]l + [Rw]) ∑ l)0
(22)
Finally, the mass balance for the two monomers is given by
dMj dt
2
) υMj )
w w (kpij + ktrij )pwi[oligomer][Mj]wVaq ∑ i)1
2
(kpij + ktrij)pi[Mj]p∫r)r ∑ i)1
rmax nuc
n j (r,t) F(r,t) dr (23)
jcr-1 where [oligomer] ) ∑n)0 [Pw]n.
4. Experimental Study 4.1. Experimental System. A schematic diagram of the VAc/BuA emulsion copolymerization system is given in Figure 1. The reactor has a capacity of 3 L, and the jacket water temperature is regulated as it flows through the jacket of the reactor for temperature control. An overhead condenser is attached at the top of the reactor, and the nitrogen is injected to maintain an inert environment in the reactor. The reaction mixture is circulated using a diaphragm pump at a fixed rate of 50 mL/min through the circulation line, in which an online densitometer is installed. To measure PSD, a portion of the reactants in the
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7231
circulation line are diverted to the capillary hydrodynamic fractionator (CHDF; Matec Applied Sciences) via a sample line into which the eluant is added for dilution. These properties, along with the reaction temperature measured by the resistance temperature detector, weights of reactants measured by the load cells, and the flow rates determined by the flowmeter, are transmitted by the distributed control system (Honeywell PlantScape system). The main computer communicates with the control computer for data acquisition and calculation of the optimal control inputs. Data in the main computer are transmitted to the control computer using Microsoft Excel Data Exchange, and Matlab (The MathWorks, Inc.) is used to compute the solid content and total number of particles. In the same way, the optimal control inputs calculated using Matlab are transmitted to PlantScape, which sends the signals to remote setpoint pumps to manipulate the flow rates of the feeds. For a more detailed description of the experimental setup and the measurement devices, one can refer to the original publication.31 4.2. Online Measurement. A first principles based method is used to calculate the solid content using the latex density measured by the online densitometer along with the amount of feed materials measured by the load cells. The total mass of the reactants in the reactor is transformed to the total volume of the reactants using the measured density of the latex:
Vm1 + Vm2 + Vp + Vw + Vsurf ) VT ) mT/Flatex
(24)
Here, the volume of water (Vw) and the volume of the surfactants (Vsurf) are known values from the measurement, while the volume of unreacted monomer i (Vmi) and the volume of the polymer in the latex (Vp) are the variables that must be determined. The volume change of mixing (∆Vmix) at the initial time is assumed to be constant during the batch time and used to correct eq 24:
∆Vmix ) VT - Vm1 - Vm2 - Vw - Vsurf
(25)
The following mass balance is derived for the relationship between the density and solid content:
Vm1Fm1 + Vm2Fm2 + VpFp )
mTm1
+
mTm2
(26)
T ) added to the where the total mass of monomer i (mmi latex up to the current time is measured by the cell load and the ratio of the conversion of the two individual monomers is computed as follows:
(VTm2 - Vm2)Fm2 (VTm1 -Vm1)Fm1 )R MWm1 MWm2
(27)
T where Vmi is the total volume of monomer i fed to the reactor until the current time. The ratio of the reaction rate of monomer 1 to that of monomer 2 (R) is calculated using the concentration of monomer i in the particle phase ([Mi]p) and the propagation rate constant (kpij) for a polymer of type i with monomer j:
R)
kp11p1[M1]p + kp21p2[M1]p kp12p1[M2]p + kp22p2[M2]p
(28)
The concentration of monomers in the particle phase (required in eq 26) is calculated from the following partitioning equation:
[Mi]pVsp +
[Mi]p [Mi]p VmiFmi Vw + Vd ) Kpi Kdi MWmi
i ) 1 and 2 (29)
where the swollen volume of the particles is calculated from the equation
Vp
Vsp )
2
1-
(30)
[Mi]pMWmi/Fmi ∑ i)1
and the volume of monomer droplets (Vd) can be calculated using the monomer concentration in the particle phases and the partition coefficients. It is assumed in this study that monomer droplets are absent in the reactor because the flow rate of monomer is very low. Consequently, eqs 24, 26, 27, and 29 are solved to calculate the five unknown variables (Vmi, Vp, [Mi]p, i ) 1 and 2), which are used to derive the solid content of the reactants. Weight-averaged PSD and the total number of particles are measured online by the CHDF. Because the CHDF measures the relative weight-averaged PSD with the largest peak assigned a value of 100, the actual distribution should be calculated through algebraic manipulation:
W(r) ) N(r)/
∑j Nj
(31)
where W(r) and N(r) represent the actual weightaveraged PSD and the relative distribution measured by the CHDF, respectively, while Nj denotes the relative distribution at discrete points of size rj. On the basis of the weight-averaged distribution, the particle density function [F(r)] is retrieved as follows:
Fi )
Vp W i NA υi
(32)
b,i b,i+1 F(r) dr] and Wi [)∫rrb,i W(r) dr] denote Here, Fi [)∫rrb,i-1 the sampled versions of the particle density function and the weight-averaged PSD, respectively, when rb,i is the upper boundary of finite element i (ri is the representative size for finite element i). NA is Avogadro’s number, and Vp is the volume of the polymer calculated in the previous density-solid content correlation equations. The volume of a particle at the discrete point of size ri is calculated as υi ) 4πri3/3, and Wi is the weightaveraged distribution at the discrete point of size ri. From the particle density function, the total number of particles can be computed as
N p ) NA
∑i Fi
(33)
Figure 2 shows the time evolution of the total number of particles and the solid content for the case of an openloop nominal input trajectory. A small deviation of the simulation results from the online measurements of the total number of particles is observed because of model mismatch coupled with a precision problem in the
7232
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
Figure 3. Open-loop nominal input trajectories. Table 1. Values of Input Levels Used in the Generation of a Database [mol/s]
Figure 2. Time evolution of (a) the total number of particles and (b) the solid content.
measurement device. Meanwhile, online measurements of the solid content compared with the simulation results (along with the gravimetric method) show negligible errors for the density-solid content calculation. In summary, the experimental results validate the accuracy of the constitutive equations with the online measurements and properties. 5. Results and Discussion In this study, the full PSD at the final batch time in a VAc/BuA emulsion copolymerization reactor is controlled using the PLS model-based predictive controller manipulating the profiles of VAc and surfactant. To reduce the dimension of the output variables (grid points for the numerical solution of the weight-averaged PSD), the Gaussian distribution
W(r,p) ) k1e-(r-a1) /σ1 + k2e-(r-a2) /σ2 2
2
(34)
is used to parametrize the distribution in such a way that the output variables used in the Y block are parameters p ) {ai, ki, σi|i ) 1, 2} for the two Gaussian distributions, that is, mean radii, amplitudes, and variances for primary and secondary modes, respectively. The total number of particles (Np) and solid content (Sc) are included in the X block to augment the information about the current state. It should be noted that, because the measurement delay for Np is 13 min (CHDF), the delayed measurement of the total number of particles is used. The inputs are the feed flow rate of
level
VAc feed flow rate (u1)
surfactant feed flow rate (u2)
1 2 3 4 5 6 7 8
0 1.495 × 10-4 2.990 × 10-4 4.485 × 10-4 5.980 × 10-4 7.475 × 10-4 8.970 × 10-4 1.047 × 10-3
0 7.10 × 10-4 1.42 × 10-3 2.13 × 10-3 2.84 × 10-3 3.55 × 10-3 4.29 × 10-3 4.97 × 10-3
VAc and the flow rate of the surfactant feed, and these inputs are also parametrized with an appropriate resolution. In the simplest manner, the input trajectory for each interval (in this study, 13 min for the measurement of Np) is chosen to be constant, which results in a single parameter for each interval. To address the nonlinearity of the system, squared inputs are included in the input data block and deviation terms from the nominal inputs (cf. Figure 3) and corresponding outputs are used because the order of magnitude changes abruptly in a batch/semibatch reactor. As a result, the X block consists of two measurements, parameters for past input and parameters for future input:
xT(t) ) [∆Np(t - td), ∆Sc(t), ∆upastT, ∆ufutureT, ∆upast,sqT] (35) Here, ∆upast,sqT ) [∆upast,12, ..., ∆upast,k2] and td is the time delay. Np and Sc represent the total number of particles and the solid content, respectively. The reactor is initially charged with deionized 1 L of water, 52 g of VAc, 0.1 g of ferrous ammonium sulfate, and 1.12 g of sodium benzoate. The reaction temperature is maintained at 67.5 °C. In the generation of a representative database for PLS model development, several methods of developing a representative ensemble are reported in the literature.32-34 In this study, the data set is generated by imposing pseudo-eight-level inputs (cf. Table 1) on the mechanistic model for 40 batches, and the resulting weight-averaged PSDs at the final time are presented in Figure 4. The setpoint is determined by imposing the nominal input trajectory shown in Figure 3 on a mechanistic model. The sampling of Np starts at 3 min, and the measurement delay is 13 min; hence, the sample time for control
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7233
Figure 4. Weight-averaged PSD at the final time generated by pseudo random eight-level signals (40 batches). Figure 6. Input trajectories calculated by a PLS model-based predictive controller and nominal input trajectories.
position of each peak more than the other factors in Gaussian distributions. Values of 0.2 and 0.4 are used in the element of the weight matrix for the current control input [∆u1(t) and ∆u2(t)], while 0.1 and 0.2 are used for future inputs. The number of latent variables (LV) is chosen as 15 to cover 91.15% of the variance in the input block and 94.67% of the variance in the output block. In general, latent variables accounting for 80% variance are sufficient for the prediction using the PLS model. However, LVs leading to more than 90% variance are used in this study to compensate for the measurement noise and highly nonlinear behavior of the process and thus to improve the predictive ability of the controller. The constraints on the input magnitude are specified as Figure 5. Weight-averaged PSD at the final time controlled by a PLS model-based predictive controller.
starts at 16 min and increases by 13 min for each subsequent sample. Therefore, the database is generated every 13 min starting from 16 min with the structure given in eq 35, and PLS is applied to the database. As shown in Figure 4, the polymerization under the conditions in this study produces a bimodal distribution: one peak by micellar nucleation (secondary peak on the left-hand side) and the other by homogeneous nucleation (primary peak on the right-hand side). The primary peaks have separate groups in which the peaks with mean values larger than 450 nm are produced when the surfactant feed is not injected during the first 16 min. Because the amount of free surfactant in the reactor is insufficient for micellar nucleation until later, micellar nucleation is delayed, and thus the primary peak by homogeneous nucleation grows larger than the peak produced under the conditions when surfactant is fed before 16 min. Figure 5 presents the weight-averaged PSD controlled by the PLS model-based predictive controller for the nominal case. The MATLAB PLS Toolbox 2.0 (Eigenvector, Inc.) is used to carry out the PLS analysis. The weight matrix for the control error is specified as Λy ) diag([1.0, 1.1, 1.0, 1.0, 1.1, 1.0]). Larger values are assigned to the elements related to the mean values of each Gaussian distribution to focus on the control of the
0 e feed flow rate of VAc [mol/s] e 8.5 × 10-4 (36) 0 e feed flow rate of surfactant [mol/s] e 4.0 × 10-3 (37) The constraints on the rate of input change are not considered because the remote setpoint pump in the experimental system changes the flow rate within 3 s, which is sufficiently short compared to the sample time. As shown in Figure 5, the primary peak in the region of large particle size is not in good agreement with its setpoint and the secondary peak in the region of small particle size is very small. This result can be explained by the database used in this case. Because the database is entirely dependent upon the mechanistic (first principles) model, the controller cannot calculate the optimal trajectory that compensates for the measurement noise and the model mismatch (cf. Figure 6). In addition, the characteristics of the primary peak reveal additional explanations. Because the primary peak is produced by homogeneous nucleation dominant in the beginning stage of the reaction (earlier than 10 min in this system29), the control actions after 16 min cannot completely eliminate the effect of disturbances (model mismatch) on this peak. To improve the utility of the database, it is updated using the experimental results under five different conditions. Figure 7 presents the weight-averaged PSD at the final time produced by conducting experiments
7234
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
Figure 7. Weight-averaged PSD at the final time generated in the experimental study (five batches).
Figure 9. Weight-averaged PSD at the final time controlled by a PLS model-based predictive controller with an updated database.
Figure 8. Input trajectories calculated by a PLS model-based predictive controller and nominal input trajectories.
Figure 10. Input trajectories calculated by a PLS model-based predictive controller with an updated database and nominal input trajectories.
using the input trajectories in Figure 8. The inputs between 0 and 16 min are specified as the nominal input values because no control action is taken in this study, and the other inputs are determined so that their values are around the nominal input values. It is noted that the setpoint is placed within the region of distributions in Figure 7. If the update of the database is introduced as the batch is repeated, the algorithm effectively introduces a batch-to-batch update.13,20 However, because the number of training batches produced by simulation is larger than the number of experimental sets, it is assumed that an improvement in the performance would be observed after five batches. Therefore, experiments are conducted in an open-loop manner under different conditions, and the experimental data are included in the existing database. Figure 9 shows the controlled distribution at the final time when the same tuning variables (the number of latent variables and weighting matrices) as those in the previous case are used. The primary peak still shows a deviation, but the secondary peak has a distribution that is closer to the setpoint. Because of possible failure in the dilution process, the particles may coagulate with each other and, thus, a trimodal distribution is observed. The bulge observed on the right-hand side of the
secondary peak also shows the influence of coagulation. Although the updated input trajectories (Figure 10) are not significantly different from those of the nominal case, the high sensitivity of the micellar nucleation in the early stage leads to a significantly different controlled response. The flow rate of the surfactant feed at 16 min is higher than the nominal value, and considering that micellar nucleation occurs after about 10 min and is sensitive to the amount of free surfactant in an earlier reaction time,29 the increase of the surfactant feed flow rate between 16 and 29 min compensates for the model mismatch. The flow rate of VAc at 16 min also explains the improved results in the secondary peak. Because the amount of monomer is related to the growth of the particle size and reduces micellar nucleation, the improved controller decreases the flow rate of VAc feed, whereas the controller based on the original simulated database increases it. Figures 11 and 12 show the weight-averaged PSD controlled by the PLS model-based predictive controller and corresponding input trajectories, respectively, when strong disturbances are introduced. The initial charge of VAc in the reactor is decreased by 50% and the concentration of the surfactant feed solution is increased
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7235 Table 2. Parameters for Two Gaussian Distributions Used in the Estimation of the Weight-Averaged PSD at the Final Time parameter k1 a1 σ1 k2 a2 σ2
unit nm nm2 nm nm2
setpoint
controlled distribution (% error)
uncontrolled distribution (% error)
9.3238 × 10-3 3.5822 × 101 3.3683 × 102 1.7063 × 10-2 1.9719 × 102 3.1265 × 103
5.3454 × 10-3 (-42.67) 2.0306 × 101 (-43.32) 8.6350 × 101 (-74.36) 2.7762 × 10-2 (62.70) 1.7139 × 102 (-13.09) 1.6297 × 103 (-47.87)
2.0332 × 10-2 (118.07) 6.6200 × 101 (84.80) 1.1794 × 103 (250.14) 6.1498 × 10-3 (-63.96) 2.5368 × 102 (28.64) 2.7533 × 103 (-11.94)
Figure 11. Weight-averaged PSD at the final time controlled by a PLS model-based predictive controller with an updated database in the presence of disturbances.
surfactant in the reactor, and the controller stops feeding the surfactant at 29 min because the disturbance in the surfactant feed is not completely eliminated between 16 and 29 min. Because the input for the surfactant feed is too small, the feed flow rate of VAc at 29 min is decreased to maintain micellar nucleation at a desired rate. Inputs for the surfactant feed after 42 min show values larger than the nominal values to compensate for aggressive input between 29 and 42 min. However, the secondary peak in controlled distribution is smaller than the desired peak, which means that the effect of aggressive input is not entirely canceled. The result shows that the effect of disturbance on the initial charge of VAc is not rejected because the initial charge error effects homogeneous nucleation at the beginning stage, and thus the mean value of the primary peak is smaller than the desired value. The performance of the controller is analyzed in a quantitative manner. Both the controlled and uncontrolled distributions are estimated using two Gaussian distributions (cf. eq 34), and the values, along with the percent relative error [(p - psp)/psp] × 100 where p ) {ki, ai, σi|i ) 1, 2}), are listed in Table 2. It is observed that mean values (a1 and a2 are mean radii [nm]) of the controlled distribution are much closer to those of the target distribution, and the amplitude and variance of the controlled secondary peak (k1 and σ1) reveal that the controller effectively eliminated the disturbances. 6. Summary
Figure 12. Input trajectories calculated by a PLS model-based predictive controller with an updated database and nominal input trajectories in the presence of disturbances.
by 26.34%. Because it is expected that the controller calculates aggressive control inputs to reject strong disturbances imposed on the system, a more robust tuning of the input weights results in 0.9 for the current input terms. The effect of disturbances is presented in Figure 11. The decrease of the initial charge of VAc causes the decrease of the growth of the primary peak, while the increase of the surfactant feed concentration produces a large secondary peak in the distribution. Although larger values were used for the control input weights than those in the previous case, the control input trajectories still show an aggressive response. The flow rate of the surfactant feed at 16 min is less than the nominal value to reduce the excess amount of
The PSD at the final time in a semibatch VAc/BuA emulsion copolymerization reactor is controlled using a PLS model-based predictive controller. The data set used to regress the statistical model is generated from a mechanistic model. The relationship between the input and output blocks is analyzed using the PLS method, and prediction models are obtained at every sample time. The prediction model at each sample time is used to update the optimal trajectory, and the strategy has been applied to an experimental system. The total number of particles with a 13 min measurement delay and solid content are used to infer the current state of the reactor, and the future input trajectories are calculated using a predictive control algorithm. Experimental results with the simulated data set show a model mismatch problem and, thus, experimental results are included in the simulated database to improve the data set. The updated database contains information about the actual system and allows the controller to calculate the optimal input trajectory to compensate for the model mismatch, and it is expected that the performance of the controller will be improved if the database is updated in a batchwise manner. However, the results show the limitation of the proposed control strategy. The homogeneous nucleation takes place dominantly in the early stage of the reaction
7236
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
(before 10 min in the present study), while the first control action is taken at 16 min. Therefore, the primary peak is hard to control, and it is recommended that, when this controller is applied, one should be cautious with respect to errors or disturbances in the early stage of the reaction because of the limited availability of the manipulation. The validity of proposed control strategy is corroborated in the presence of strong disturbances. The disturbance in the surfactant feed was eliminated satisfactorily, while the initial charge error could not be rejected because of the characteristics of the system. Consequently, the proposed control strategy has been proven to be effective for the regulation of PSD in a semibatch emulsion copolymerization reactor. Acknowledgment This work was partially supported by the Postdoctoral Fellowship Program of the Korea Science and Engineering Foundation (KOSEF). Nomenclature Cmicelle ) concentration of micelles D0 ) diffusion coefficient l ) entry rate constant for oligomers of type i ei,micelle H ) heaviside function, which is unity when the argument is nonnegative and zero otherwise kd1 ) kinetic constants for the oxidation step kw pav ) pseudohomogeneous rate constant for propagation 2 2 w in the aqueous phase (kw pav ) ∑i)1 ∑j)1 kpijpwi[Mj]w) kpij ) rate constant for propagation of a polymer of type i with monomer j w kpij/kpij ) rate constants in the particle/aqueous phases for propagation for a polymer of type i with monomer j kri ) rate constant for reaction between initiator and monomer i kw tav ) pseudohomogeneous rate constant for termination 2 2 w in the aqueous phase (kw tav ) ∑i)1 ∑j)1 ktijpwipwj) w ktrij/ktrij ) rate constants in the particle/aqueous phases for chain transfer to monomer for a polymer of type i with monomer j w ktrij/ktrij ) rate constants in the particle/aqueous phases for chain transfer to monomer for a polymer of type i with monomer j [Mi]w ) concentration of monomer i in the aqueous phase [Mj]p ) concentration of monomer j in the particles MWj ) molecular weight of monomer j T mmi ) total mass of monomer i added to the latex mT ) total mass of the latex n j (r,t) ) average number of active radicals in a particle of size r at time t NA ) Avogadro’s number pi ) probability that a radical is of type i in the particles [Pw]l ) concentration of oligomer with chain length l in the aqueous phase pwi ) probability of a radical of type i in the aqueous phase rcutoff ) cutoff size below which the particles are prone to coagulate rI ) stoichiometric ratio between the oxidizer and reducer rupper ) maximum size of the particles that could result by the coagulation of smaller particles Vaq ) volume of the aqueous phase vi ) molar feed rate of component i Vmi ) volume of monomer i T ) total volume of monomer i fed to the reactor Vmi Vp ) volume of the polymer in the latex Vsp ) swollen volume of the particles
Vsurf ) volume of the surfactants VT ) total volume of the latex Vw ) volume of water WS ) Fuch’s stability ratio [Yr1] ) concentration of the catalyst in the reduced form Greek Letters β ) coagulation kernel Flatex ) density of the latex Fmi ) density of monomer i Fp ) density of polymer
Literature Cited (1) Ramkrishna, D. Population Balances; Academic Press: San Diego, 2000. (2) Park, M.-J.; Doyle, F. J., III. Nonlinear Model Order Reduction and Control of Particle Size Distribution in a Semibatch Vinyl Acetate/Butyl Acrylate Emulsion Copolymerization Reactor. Korean J. Chem. Eng. 2004, 21, 168-176. (3) Kalani, A.; Christofides, P. D. Simulation, Estimation and Control of Size Distribution in Aerosol Processes with Simultaneous Reaction, Nucleation, Condensation and Coagulation. Comput. Chem. Eng. 2002, 26, 1153-1169. (4) Zeaiter, J.; Romagnoli, J. A.; Barton, G. W.; Gomes, V. G.; Hawkett, B. S.; Gilbert, R. G. Operation of Semi-batch Emulsion Polymerisation Reactors: Modeling, Validation and Effect of Operating Conditions. Chem. Eng. Sci. 2002, 57, 2955-2969. (5) Valappil, J.; Georgakis, C. Nonlinear Model Predictive Control of End-Use Properties in Batch Reactors. AIChE J. 2002, 48, 2006-2021. (6) Sharaf, M. A.; Illman, D. L.; Kowalski, B. R. Chemometrics; John Wiley & Sons: New York, 1986. (7) Geladi, P.; Kowalski, B. R. Partial Least-Squares Regression: A Tutorial. Anal. Chim. Acta 1986, 185, 1-17. (8) Dayal, B. S.; MacGregor, J. F. Improved PLS Algorithms. J. Chemom. 1997, 11, 73-85. (9) Lindgren, F.; Geladi, P.; Wold, S. The Kernel Algorithm for PLS. J. Chemom. 1993, 7, 45-59. (10) Mejdell, T.; Skogestad, S. Estimation of Distillation Compositions from Multiple Temperature Measurements Using PartialLeast-Squares Regression. Ind. Eng. Chem. Res. 1991, 30, 25432555. (11) Durand, J.-F. Local Polynomial Additive Regression through PLS and Splines: PLSS. Chemom. Intell. Lab. Syst. 2001, 58, 235246. (12) Baffi, G.; Martin, E. B.; Morris, A. J. Non-linear Projection to Latent Structures Revisited: The Quadratic PLS Algorithm. Comput. Chem. Eng. 1999, 23, 395-411. (13) Flores-Cerillo, J.; MacGregor, J. F. Within-Batch and Batch-to-Batch Inferential-Adaptive Control of Semibatch Reactors: A Partial Least Squares Approach. Ind. Eng. Chem. Res. 2003, 42, 3334-3345. (14) Prautzsch, H.; Boehm W.; Paluszny, M. B’ezier and Bspline techniques; Springer: New York, 2002. (15) Wold, S.; Antti, H.; Lindgren, F.; O ¨ hman, J. Orthogonal Signal Correction of Near-infrared Spectra. Chemom. Intell. Lab. Syst. 1998, 44, 175-185. (16) Trygg, J.; Wold, S. Orthogonal Projections to Latent Structures (O-PLS). J. Chemom. 2002, 16, 119-128. (17) Wold, S.; Kettaneh-Wold, N.; Skagerberg, B. Non-linear PLS Modelling. Chemom. Intell. Lab. Syst. 1989, 7, 53-65. (18) Holcomb, T. R.; Morari, M. PLS/Neural Networks. Comput. Chem. Eng. 1992, 16, 393-411. (19) Elizalde, O.; Vicente, M.; Leiza, J. R.; Asua, J. M. Control of the Adhesive Properties of n-Butyl Acrylate/Styrene Latexes. Polym. React. Eng. 2002, 10, 265-283. (20) Doyle, F. J., III; Harrison, C. A.; Crowley, T. J. Hybrid Model-Based Approach to Batch-to-Batch Control of Particle Size Distribution in Emulsion Polymerization. Comput. Chem. Eng. 2003, 27, 1153-1163. (21) Kesavan, P.; Lee, J. H.; Saucedo, V.; Krishnagopalan, G. A. Partial Least Squares (PLS) based Monitoring and Control of Batch Digesters. J. Process Control 2000, 10, 229-236. (22) Belsley, D.; Kuh, E.; Welsch, R. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity; Wiley: New York, 1980.
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7237 (23) Flores-Cerillo, J. Quality Control for Batch Processes Using Multivariate Latent Variable Methods. Ph.D. Thesis, McMaster University, Ontario, Canada, 2003. (24) Nelson, P. R. C.; Taylor, P. A.; MacGregor, J. F. Missing Data Methods in PCA and PLS: Score Calculations with Incomplete Observations. Chemom. Intell. Lab. Syst. 1996, 35, 45-65. (25) Arteaga, F.; Ferrer, A. Dealing With Missing Data in MSPC: Several Methods, Different Interpretations, Some Examples. J. Chemom. 2002, 16, 404-418. (26) Johansen, T. A.; Foss, B. A. Operating Regime Based Process Modeling and Identification. Comput. Chem. Eng. 1997, 21, 159-176. (27) Stenman, A. Model on Demand: Algorithms, Analysis and Applications. Ph.D. Thesis, Linko¨ping University, Linko¨ping, Sweden, 1999. (28) Hur, S.-M.; Park, M.-J.; Rhee, H.-K. Design and Application of Model-on-Demand Predictive Controller to a Semibatch Copolymerization Reactor. Ind. Eng. Chem. Res. 2003, 42, 847859. (29) Immanuel, C. D.; Cordeiro, C. F.; Sundaram, S. S.; Doyle, F. J., III. Population Balance PSD Model for Emulsion Polymerization with Steric Stabilizers. AIChE J. 2003, 49, 1392-1404. (30) Immanuel, C. D.; Doyle, F. J., III. ComputationallyEfficient Solution of Population Balance Models Incorporating
Nucleation, Growth and Coagulation: Application to Emulsion Polymerization. Chem. Eng. Sci. 2003, 58, 3681-3698. (31) Immanuel, C. D.; Crowley, T. J.; Meadows, E. S.; Cordeiro, C. F.; Doyle, F. J., III. Evolution of Multimodal Particle Size Distribution in Vinyl Acetate/Butyl Acrylate Emulsion Copolymerizations. J. Polym. Sci., Part A: Polym. Chem. 2003, 41, 22322249. (32) Bangia, A. K.; Batcho, P. F.; Kevrekidis, I. G.; Karniadakis, G. E. Unsteady 2-D Flows in Complex Geometries: Comparative Bifurcation Studies with Global Eigenfunction Expansions. SIAM J. Sci. Comput. 1997, 18, 775-805. (33) Graham, M.; Kevrekidis, I. G. Alternative Approaches to Karhunen-Loeve Decomposition for Model Reduction and Data Analysis. Comput. Chem. Eng. 1996, 20, 495-506. (34) Loffler, H. P.; Marquardt, W. On the Order Reduction of Nonlinear Differential-Algebraic Process Models. Proc. Am. Control Conf. 1992, 1546-1551.
Received for review December 14, 2003 Revised manuscript received March 15, 2004 Accepted April 2, 2004 IE034304G
