Control of batch product quality by trajectory manipulation using latent variable models

Journal of Process Control 14 (2004) 539–553 www.elsevier.com/locate/jprocont

Control of batch product quality by trajectory manipulation using latent variable models Jesus Flores-Cerrillo, John F. MacGregor

*

Department of Chemical Engineering, McMaster University, Hamilton, Ontario, Canada L8S4L7 Received 25 June 2003; received in revised form 22 September 2003; accepted 22 September 2003

Abstract A novel inferential strategy for controlling end-product quality properties by adjusting the complete trajectories of the manipulated variables is presented. Control through complete trajectory manipulation using empirical models is possible by controlling the process in the reduce space (scores) of a latent variable model rather than in the real space of the manipulated variables. Model inversion and trajectory reconstruction is achieved by exploiting the correlation structure in the manipulated variable trajectories captured by a partial least squares model. The approach is illustrated with a condensation polymerisation example for the production of nylon and with data gathered from an industrial emulsion polymerisation process. The data requirements for building the model are shown to be modest.
2003 Elsevier Ltd. All rights reserved. Keywords: Product quality; Partial least squares; Reduced space control

1. Introduction Batch/semi-batch processes are commonly used because their flexibility to manage many different grades and types of products. In these processes, it is necessary to achieve tight final quality specifications. However, this is not easily achieved because batch operations suffer from constant changes in raw material properties, variations in start-up initialisation, and in operating conditions, all of which introduce disturbances in the final product quality. Moreover, compensating for these disturbances is difficult due to the non-linear behaviour of the chemical reactors and to the fact that robust online sensors for monitoring quality variables are rarely available. Control of product quality usually requires the online adjustment of several manipulated variable trajectories (MVTs) such as the pressure and temperature trajectories. Several approaches based on detailed theoretical models have been presented. These can generally be divided into two groups, the first based on

*

Corresponding author. Tel.: +1-905-525-9140x24951; fax: +1-905521-1350. E-mail address: [email protected] (J.F. MacGregor). 0959-1524/$ - see front matter
2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jprocont.2003.09.008

non-linear differential geometric control, and the second based on on-line optimization. The differential geometric approaches [1–3] use the non-linear model to perform a feedback transformation that linearizes the system and then linear control theory can be applied. Examples in the literature include the control of final latex properties such as instantaneous copolymer composition, conversion and weight average molecular weight common in the emulsion polymerisation of styrene–butadiene [1], and the control of copolymer composition and weight average molecular weight for the free radical polymerisation of vinyl acetate/methyl methacrylate reaction [2]. In on-line optimization, optimal trajectories are periodically recomputed at various instances throughout the batch to optimize some final quality and/or performance measure. Some examples include Crowley and Choi [4] for the on-line control of molecular weight distribution and conversion on the free radical polymerisation of methyl methacrylate, and Ruppen et al. [5] for on-line batch time minimization and conversion control in an experimental set-up. In both approaches control action was obtained using sequential quadratic programming methods at several time intervals. In spite of the significant literature addressing the trajectory control of batch processes, many of these

540

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

Nomenclature number of principal components residual matrix number of on-line measurements for the jth variable F residual matrix g number of off-line analysis for the sth variable K number of batches l number of trajectories for the on-line variables M number of quality properties n number of trajectories for the manipulated variables PT loading matrix pT loading vector Q1 weighting matrix in the controlled scores Q2 score suppression movement matrix QT projection matrix from PLS r number of the off-line variables s2 variance of a score T score matrix t score vector ^tpresent vector of estimated scores uc vector of manipulated variables trajectories uc;future vector of future control actions (hi 6 h 6 hf ) uc;implemented vector of implemented control actions (0 6 h 6 hi
1 ) w number of segments for the mth manipulated variable W projection matrix A E f

strategies are difficult to implement because they are computationally intensive and/or require substantial model knowledge. Recently, Bonvin and co-workers [6–8], recognizing that the use of detailed theoretical models for the control and optimization of batch processes is unrealistic in industry, introduce a strategy in which the optimal structure of the parameterised inputs is determined using, for example an approximate model and then measurements (off-line and/or off-line) are employed to refine (update) them. Empirical modelling, on the other hand, has the advantage of ease in model building. Yabuki and MacGregor [9,10], and Flores-Cerrillo and MacGregor [11,12] among others used empirical models for the control of product quality-properties, but in these approaches the control action was restricted to only a few movements in the manipulated variables (injection of additional reactants) because, in these cases, these few adjustments were enough to reject the disturbances and to achieve the desired end-qualities. However, if the operation calls for adjustments to MVTs through most of the duration of the process, another approach needs

regressor vector that includes on-line and offline measurements, and control actions X unfolded regressor matrix of process trajectories (MVTs and measurements) v three-dimensional array vector of total measurements (on-line and xm off-line) xm;future vector of unmeasured variables at time hi (hiþ1 6 h 6 hf ) xm;measured vector of measured variables at time hi (0 6 h 6 hi ) xoff vector of off-line measurements xon vector of on-line trajectory measurements Y matrix of quality properties y vector of quality variables ^y vector of estimated quality variables x

Greek symbols k weighting factor h decision times d de-tuning factor a proportionality vector b coefficients for the PLS inner relation Index a i j,s,m k f

latent variable index time index variable index batch index final batch time

to be taken. The approach often used in these cases is to segment the MVTs into a small number of intervals (e.g. 5–10) and force the behaviour of the MVTs over the duration of each interval to follow a zero or first order hold. Control is then accomplished by manipulating the slope or the level (stair-case parameterisation) at the start of each interval (decision points). Studies involving this type of parameterisation can be found in [13,14] among others. However, in many batch processes such a staircase parameterisation of the MVTs, just for convenience of the control engineers, may not be acceptable. The operation of the batch may require, or historically be based on, smooth MVTs, and converting them to stair-case approximations might represent a radical departure from normal practice, with the implication that control schemes based on them will never be implemented. Moreover, model inversion in the control algorithm would be usually difficult with this approach because a large number of highly correlated control actions need to be determined at every decision point. A solution to this problem comes from recognizing that within the range of normal process operation all the

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

process variable trajectories (both MVTs and measured variables) are very highly correlated with one another, both contemporaneously (i.e. at the same time period) and temporally (over the time history of the batch). This implies that their behaviour can be represented in a much lower dimensional space using latent variable models based on principal component analysis (PCA) or partial least squares (PLS). This concept has been powerfully exploited for the analysis and monitoring of batch processes [15–17] where the entire time histories of all the process and MVTs can usually be summarized by only a few (2 or 3) latent variables. Therefore, in this paper we show that by projecting all the process variable trajectory data into low dimensional latent variable spaces, all control decisions can be performed on the latent variables, and the entire MVTs for the remainder of the batch then reconstructed from the latent variable models. In this reduced dimensional space, the data requirements for modelling and for model parameter estimation are much less demanding, the control computation is easier, and the computed MVTs are smooth and consistent with past operation of the process. In spite of these inherent advantages in controlling the MVTs of batch processes in a latent variable space, no literature has yet addressed this issue. Reduced dimension controllers for continuous processes (a binary distillation column simulator and the Tennessee Eastman process) based on PCA have been proposed [18–20] which express the control objective in the score space of a PCA model, but the dimension of the manipulated variable space is still small since no trajectories need to be computed. The purpose of this paper is to introduce an inferential control strategy that allows a much finer characterisation and smoother reconstruction of optimal MVTs than those obtained using staircase parameteri-

541

sation, and one that reduces the complexity and number of identification experiments needed for model building. These objectives are made possible by formulating the control strategy in the reduced dimensional space of a latent variable model, and then inverting the model to obtain the solution for the MVTs. The outline of the paper is as follows: in Section 2 the methodology is introduced; in Section 3, the control approach is illustrated with a condensation polymerisation case study for the production of nylon and preliminary results are shown for an industrial emulsion polymerisation process.

2. Control methodology 2.1. Model building The proposed methodology uses historical databases and a few complementary identification experiments for model building. The empirical model is obtained using PLS. However, other projection methods such as principal component regression may also be applied. The database from which the PLS model is identified is shown in Fig. 1. It consists of a (K  M) response matrix Y and an originally three-dimensional array v, which after unfolding [17,22] would yield a (K  N ) regressor matrix X where K is the number of batches. Each row vector of Y denoted as yT , contains M quality properties measured at the end of each batch. Each row vector of X, denoted as xT , is composed of:
 xT ¼ xTon xToff uTc where xTon ¼ ½ xTon;1 xTon;2    xTon;l  is a vector of the trajectories of l on-line process variables such as temperature and pressure obtained from on-line sensors;

Fig. 1. Unfolding of database for model building.

542

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

xToff ¼ ½ xToff;1 xToff;2    xToff;r  is the set of any off-line measurements collected occasionally on r variables during the batch, and uTc ¼ ½ uTc;1 uTc;2    uTc;n  is a vector of the trajectories of n manipulated variables. As can be seen in Fig. 1, xTon;j ¼ ½xon;1 ; . . . xon;f j and xToff;s ¼ ½xoff;1 ; . . . ; xoff;g s denotes, respectively, the row vector of observations obtained from on-line measurements on the jth variable, and from off-line measurements on the sth variable over the course of the batch, while uTc;m ¼ ½ uc;1 ;    uc;w m denotes the trajectory of the mth manipulated variable (MV). Here, f , g and w are, respectively, the number of on-line measurements, off-line analysis and MV segments for the corresponding variable in each category. Therefore the regressor matrix is of dimension (K  N ) where N ¼ fl þ gr þ wn. In the following text xTon and xToff are combined into a single row vector xTm ¼ ½ xTon xToff , and then xT ¼ ½ xTm uTc . Full MVTs are obtained through trajectory segmentation as illustrated in Fig. 2. The MVTs are segmented into a (possibly) large number of intervals ðwÞ and control decision points (hi ; i ¼ 1; 2; . . .) are selected. At each decision point (hi ), final properties (y) are predicted and the adjustments to the remaining MVTs (after this decision point) are computed if the predicted final properties are not within desired specifications. Notice that the segment size is not necessarily uniform and that decisions points may be chosen arbitrarily but are assumed to be the same for each batch. (The decision points will usually be selected using prior process knowledge.) In the limit, control action can be taken at every segment (i.e. every segment would represent a decision point), but this is almost never necessary, as a very small number is usually adequate. The fineness of the trajectory segmentation will largely depend on how fine the shape of the trajectories needs to be reconstructed. The control methodology presented in the paper is essentially independent of this. The data-set used for model building consists of representative operating data from past batches in order to capture information on most of the disturbances and operating policies normally encountered in the batch

Fig. 2. Fine segmentation of MVTs and decision points.

operation. In addition, data in which some changes in the MVTs are performed at each decision point are required in order to establish causal relationship between these MVT changes and the other measured process variable trajectories and the final product qualities. Rebuilding the model by adding new batch data collected after implementing the control scheme can also be done in order to further improve the causal relationship and expand the information on the effect of disturbances on the trajectories. The data requirements are further discussed in the examples. Linear PLS regression is then performed by projecting the scaled (unit variance) data (expressed as deviations from their nominal conditions) onto lower dimensional subspaces: X ¼ TPT þ E Y ¼ TQT þ F

ð1Þ

where the columns of T are values of new latent variables (T ¼ XW) that capture most of the variability in the data, P and Q are the loading matrices for X and Y respectively, and E and F are residual matrices. Nonlinear PLS regression can also be used as will be shown at the end of Section 3.1. However, for simplicity, in the following discussion linear models are assumed. The control methodology used in this work consists of two stages: at predetermined decision times (hi , i ¼ 1; 2; . . .) an inferential end-quality prediction using on-line and possible off-line process measurements (xm ) and MVTs (uc ) available up to that time is performed to determine whether or not the controlled end-qualities (y) fall outside a pre-determined ‘‘no-control’’ region, and then if needed, control action is computed in the latent variable space followed by model inversion to obtain the modified MVTs for the remainder of the batch that will yield the desired final qualities. This two-stage procedure is repeated at every decision point (hi ) using all available measurements on the process variable and MVTs available up to that time. The novelty of the proposed approach is that the control and the model inversion stage is performed in the reduced dimensional space (latent variable or score space) of a PLS model rather than in the real space of the MVTs. Due to the high correlation of measurements and control actions, the true dimensionality of the process, determined by the score variable space (ta ; a ¼ 1; 2; . . . ; A) of the PLS model, is generally much smaller than the number of manipulated variable points obtained from the MVT segmentation (uc ). Therefore, the control computation performed in the reduced latent variable space (t) is much simpler than the one performed in the real space. In the following, the control methodology is described for one control decision point (hi ) during the batch. This is simply repeated at each future decision point. Notice that although the method is illustrated with an example in which the decision points are defined at fixed clock

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

times (hi ; i ¼ 1; 2; . . .), these decision points could easily be based on measured variables other than time, such as specified values of conversion or energy production. This would be an advantage on batches that do not have the same duration (due to, for example, seasonal variations in cooling capacity and varying row material properties), since the process trajectories can then be aligned using such indicator variables [21,15,22–24]. 2.2. Prediction For on-line end-quality estimation (^ y), when a new batch k is being processed, at every decision point (hi ; i ¼ 1; 2; . . .) 0 6 hi 6 hf , there exists a regressor row vector xT composed of at least the following variables:
 xT ¼ xTm uTc h i ¼ xm;measured;hTi xTm;future uc;implemented;hTi ; uTc;future ð2Þ The regressor vector x consists of: all measured variables (xm;measured ) available up to time hi (0 6 h 6 hi ); unmeasured variables (xm;future ) not available at hi , but that will be available in the future (hiþ1 6 h 6 hf ); implemented control actions uc;implemented (0 6 h 6 hi
1 ); and future control actions uc;future , (hi 6 h 6 hf ) which will be determined through the control algorithm. Note that at the model building stage, the xm;future and uc;future vectors are available for each batch. To estimate whether or not the final quality properties for a new batch will lie within an acceptable region, the prediction is performed considering uc;future ¼ uc;nominal (i.e. assuming that the remaining MVTs will be kept at their nominal conditions) using the PLS model:
 ^tTpresent ¼ xTm uTc W i h ¼ xTm;measured;hi ; xTm;future uTc;implemented;hi ; uTc;nominal W ð3Þ ^ yT ¼ ^tTpresent Q

T

ð4Þ

W and Q are projection matrices obtained from the PLS model building stage. The vector of scores, ^tpresent , for the new batch is the projection of the x vector onto the reduced dimension space of the latent variable model at time hi , and ^ y is the vector of predicted end-quality properties. From the above equations, it can be noticed that changes in batch operation detected by measurements of the process variable trajectories (xm;measured;hi ) or produced by changes in the MVTs (uc;implemented;hi ) would produce changes in the scores (^tpresent ) and therefore in the end-quality properties (i.e. changes in the endqualities can be detected through changes in the scores). From Eq. (3), it can also be noticed that in order to compute ^tpresent and ^ y, it is necessary to have an estimate of the unknown future measurements (xm;future ) from

543

(hiþ1 6 h 6 hf ). These can be imputed from the PLS model for the batch process using efficient missing data algorithms available in the literature [25,26]. Alternatively, a multi-model approach in which different PLS models are identified at every decision point can be used [14] or a recursive Kalman filter approach as shown in [14] taken. In this paper a single PLS model is used for prediction and control, and the estimation of unknown future measurements is performed using the PLS model and a missing data algorithm. Missing data imputation based on, for example, conditional expectation or expectation/maximisation (EM) have been shown to provide very powerful time-varying model predictive forecast of the remaining portions of the batch trajectories [27]. Such efficient predictions are possible because the latent variable models based on PLS (or PCA) capture the time varying covariance structure of the data over the entire batch trajectory. These predictions will be much better than those provided by fixed time series or Kalman filter models [27]. The ‘‘no-control region’’ can be determined in several ways, such as one that takes into account the uncertainty of the model for prediction [9], using product specifications, or with quality data under normal (‘‘incontrol’’) operating conditions [12]. In this work a simple control region based on product quality specifications will be used (Section 3). The issue of whether or not to use a ‘‘no-control’’ region is at the discretion of the user, and is not essential to the control methodology presented in this paper. If the quality prediction is outside the ‘‘no-control’’ region, then a control action, and model inversion to obtain the MVTs for the remainder of the batch uTc;future is needed. Obtaining the full MVTs consist of two stages: (1) computation of the adjustments required in the latent variable scores Dt, followed by (2) model inversion of the PLS model to obtain the real MVTs for the remainder of the batch. These two stages are explained in the following sections. 2.3. Score adjustment computation At every decision point (hi ), the change in the scores (Dt) needed to track the end-qualities closer to their setpoints (ysp ) can be obtained by solving the quadratic objective: min |{z}

T

ð^y
ysp Þ Q1 ð^y
ysp Þ þ DtT Q2 Dt þ kT 2

Dtðhi Þ

st

^yT ¼ ðDt þ ^tpresent ÞT QT 2 A X ðDt þ ^tpresent Þa T2 ¼ 2 sa a¼1

ð5Þ

Dtmin 6 Dt 6 Dtmax where DtT ¼ tT
^tTpresent , Q1 is a diagonal weighting matrix defining the relative importance of the variables

544

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

y’s, Q2 is a diagonal movement suppression matrix that is used as a tuning matrix to moderate the aggressiveness of the control, T 2 is the Hotelling’s statistic, s2a is the variance of the score ta , and k is a weighting factor which determines how tightly the solution is to be constrained to the region of the score space defined by past operation. Russell et al. [14] used a similar constraint on T 2 . Hard constraints in the adjustment to the scores (Dtmin 6 Dt 6 Dtmax ) are problem dependent and may or not need to be included. Soft constraints on Dt are contained in the quadratic objective function. The soft constraint on the score magnitudes through, Hotelling’s T 2 statistic, is intended to constrain the solution in the region where the model is valid. Eq. (5) is a quadratic programming problem that can be restated as: min |{z} Dtðhi Þ

1 T Dt HDt þ f T Dt 2

ð6Þ

where H ¼ Q Q1 Q þ Q2 þ Q3 Q3 ¼ diag½k=s2a  Dtmin 6 Dt 6 Dtmax

ð7Þ

In the case of no hard constraints, the solution is easily obtained as: DtT ¼
f T H
1

ð8Þ

The aim of Eq. (8) is to obtain the change in the scores (Dt) that would drive the final quality variables closer to their desired set-points (ysp ). Due to the movement suppression matrix (Q2 ) and/or k, the computed (Dt) may not drive the process all the way to their set-points. Choosing Q1 ¼ I, Q2 ¼ 0 and k ¼ 0, gives the minimum variance controller, which, at each decision point would force the predicted qualities (^ y) to be equal to their set-points (^ y ¼ ysp ) at the end of the batch: min |{z}

T

ð^ y
ysp Þ Q1 ð^ y
ysp Þ

ð11Þ

3. dimðDtÞ > dimðysp Þ This case is a common situation. Although the number of variables to be used in the control algorithm has been reduced to A latent variables, a projection from a lower to higher space is still required. In this situation Eq. (9) has an infinite number of solutions. Therefore, a natural choice is to select the Dtðhi Þ having the minimum norm-2: min |{z}

DtT Dt

st

yTsp

Dtðhi Þ

ð12Þ ¼ ðDt þ ^tpresent ÞT QT

and whose solution can be easily obtained as:
1

DtT ¼ ðyTsp
^tTpresent QT ÞðQQT Þ Q

ð13Þ


1 DtT ¼ dðyTsp
^tTpresent QT ÞðQQT Þ Q

ð14Þ

This is a simple alternative to using the quadratic term DtT Q2 Dt in the general linear quadratic control objective (5). A Dt vector is computed at every decision point (hi ). Eqs. (10), (11) and (13) are consistent with the PLS model inversion results found in [28]. Notice that in this last situation (Eq. (14)), the matrix QQT has dimension m  m (m being the number of quality properties). Therefore, in order to avoid illconditioned matrix inversion, the quality properties should not be highly correlated. This poses no problem since one can always perform a PCA on the Y quality matrix to obtain a set of orthogonal variables (s) that can be used as new controlled variables. Alternatively, if it is decided to retain an independent set of physical y variables, selective PCA [28] can be performed on the Y matrix to determine that subset of quality variables which best defines the Y space.

ð9Þ

Dtðhi Þ

st


1

A detuning factor (0 6 d 6 1) may be included for this reduced space controller in order to moderate the aggressiveness of the control moves:

T

f T ¼ ðQ^tpresent
ysp ÞT Q1 Q þ ^tTpresent Q3

DtT ¼ ðyTsp
tTpresent QT ÞQðQT QÞ

T ^ y ¼ ðDt þ ^tpresent Þ QT

2.4. Inversion of PLS model to obtain the MVTs

T

Three situations arise (for the unconstrained case) in finding a solution to (9) depending on the statistical dimensions of ysp and (Dt): 1. dimðDtÞ ¼ dimðysp Þ In this situation a unique solution exists that can be directly obtained from (9): DtT ¼ yTsp ðQT Þ


1


tTpresent

2. dimðDtÞ < dimðysp Þ In this case a least square solution is needed:

ð10Þ

Once the low dimensional (A  1) vector Dt is computed via one of the control algorithms described in the last section, it remains to reconstruct from tT ¼ DtT þ ^tTpresent , estimates for the high dimensional trajectories for the future process variables (xm;future ) and for the future manipulated variables (uc;future ) over the remainder of the batch. These future trajectories can be computed from the PLS model (1) in such a way that their covariance structure is consistent with past operation. If there were no restrictions on the trajectories, such as might be the case for a control action at h ¼ 0,

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

then the model for the X-space can be used directly to compute the x vector trajectory (xT ¼ ½ xTm uTc ) for the entire batch [28] as: xT ¼ tT PT

ð15Þ

However for control intervals at times hi > 0 the x vector trajectory (xT ¼ ½xTm;measuredð0:hi Þ uTc;implementedð0:hi Þ xTm;futureðhi :hf Þ uTc;futureðhi :hf Þ ) is composed of measured process variables (xTm;measuredð0:hi Þ ) for the interval 0 6 h < hi , and for the already implemented manipulated variables (uTc;implementedð0:hi Þ ) that must be respected when computing the trajectories for the remainder of the batch (hi 6 h < hf ). Denote xT1 ¼ ½ xTm;measuredð0:hi Þ uTc;implementedð0:hi Þ  the known trajectories over the time interval (0 : hi ) that must be respected, xT2 ¼ ½ xTm;futureðhi :hf Þ uTc;futureðhi :hf Þ  the remaining trajectories to be computed, and PT1 and PT2 their corresponding loading matrices. At times hi > 0, if x is directly reconstructed using (15) as xT ¼ tT PT then
T 
 ð16Þ x1 xT2 ¼ tT PT1 tT PT2 However, the computed tT PT1 will not be equal to the actually observed trajectories at time hi xT1 ¼ ½ xTm;measuredð0:hi Þ uTc;implementedð0:hi Þ . Therefore, simply selecting xT2 ¼ tT PT2 would not be correct as it does not account for what has actually been observed for xT1 in the first part of batch. Therefore, assume that the remaining trajectories (future manipulated variables and measurements) are: xT2 ¼ ðtT þ aT ÞPT2

ð17Þ

where aT PT2 is an adjustment to xT2 that accounts for the effects of discrepancy between tT PT1 and xT1 during the first part of the batch. (Selection of such a relationship will also ensure that the correlation structure of the PLS model is kept.) However, we still wish to achieve the computed value in score space t that will satisfy the overall PLS model. Therefore, we must have:  
T  W1 T T t ¼ x1 x2 ð18Þ ¼ xT1 W1 þ xT2 W2 W2 then xT2 W2 ¼ tT
xT1 W1

ð19Þ

Substituting xT2 ¼ ðtT þ aT ÞPT2 in (19): ðtT þ aT ÞPT2 W2 ¼ tT
xT1 W1 Therefore ðtT þ aT Þ ¼ ðtT
xT1 W1 ÞðPT2 W2 Þ


1

ð20Þ

And by substituting (20) in (17) the remaining MVTs to be implemented are obtained (hi 6 h < hf ):
1

xT2 ¼ ðtT
xT1 W1 ÞðPT2 W2 Þ PT2

ð21Þ

545

It is easy shown that this equation reduces to the relationship in (15) when hi ¼ 0 where there are no existing trajectory measurements or manipulated variables. The (A  A) matrix PT2 W2 is nearly always well conditioned, and so there is no problem with performing the inversion [29]. This inferential control algorithm is then repeated at every decision point (hi ) until completion of the batch.

3. Case studies 3.1. Case study 1. Condensation polymerisation In the batch condensation polymerisation of nylon 6,6 the end product properties are mainly affected by disturbances in the water content of the feed. In plant operation, feed water content disturbances occur because a single evaporator usually feeds several reactors [30]. The non-linear mechanistic model of nylon 6,6 batch polymerisation used in this work for data generation and model performance evaluation was developed in [30]. The complete description of the model and model parameters can be found in the original publication. This system was studied in [14,30], where several control strategies including conventional control (PID and gain scheduled PID), non-linear model based control and empirical control based on linear state-space models were evaluated. In the databased approach [14], control of the system was achieved by reactor and jacket pressure manipulation. These two manipulated variables were segmented and characterised by slope and level (stair-case parameterisation) leading to 10 control variables. A total of 7 intervals (decision points) were used. An empirical state space model was identified from 69 batches arising from an experimental design in the 10 manipulated variables. Several differences between the control strategy proposed here and the one used in [14] can be noticed, the most important being: (i) the control is computed in the reduce latent variable space rather than in the real space of the MVTs, (ii) only two decision points are needed to achieve good control; thereby simplifying the implementation and decreasing the number of identification experiments needed to build a model, and (iii) a much finer MVT reconstruction is achieved. Control objectives and trajectory segmentation: The control objective is to maintain the end-amine concentration (NH2 ) and the number average molecular weight (MWN) at their set-points to produce nylon 6,6 when the system is affected by changes in the initial water content (W). The MVTs used to control the end-qualities are the jacket and reactor pressure trajectories. These MVTs are finely segmented every 5 min starting at 35 min from the beginning of the reaction until 30 min

546

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

before the completion of the batch (total reaction time 200 min), giving trajectories defined at 40 discrete time points in the interval (35 6 h 6 170). The trajectories for the first 35 min and the last 30 min were fixed for all batches. Two control decision points at 35 and 75 min were found to be sufficient for good control for the conditions used in this example. In order to predict NH2 and MWN, on-line measurements of the reactor temperature (Tr ) and venting (v) are considered available every two minutes. Data generation: A PLS model with five latent variables (determined by cross-validation) was built from a data set consisting of 45 batches in which the initial water content (W ) was randomly varied. In 30 of the batches some movement in the MVT (at the two decision points) was performed (some of these batches would normally be available from historical data). The effect that the number batches used for identification of the PLS model has on control performance is discussed at the end of this section and in [29]. Prediction: The first step is to evaluate the performance of the PLS model prediction with different missing data algorithms at each decision point. Several missing data algorithms were tested and as an illustration some results are shown in Figs. 3 and 4. The predicted trajectories, made using the available data up to the first decision point (hi ¼ 35 min), for venting (v) and reactor temperature (Tr ) when the process is affected by a disturbance of )10% (mass) in the initial water content are shown in Fig. 3. Each predicted trajectory is obtained using: (-·-) expectation–maximisation (EM), (-h-) iterative-imputation (IMP), (  ) single component projection (SCP), and (--) projection to the plane using PLS (PTP) method [26,27]. As judged from this example and many similar simulations, all the missing data algorithms provide reasonable estimates of the trajectories, except perhaps the SCP method. In Fig. 4 predictions of the final qualities made at the first decision point (hi ¼ 35 min) using the IMP approach are shown when the initial water content randomly varies for 15 batches in the range of ±10% (mass). As can be seen in this figure, the predicted final quality properties (at h ¼ 200 min) made using the PLS model at the first decision point (hi ¼ 35 min) are in good agreement with the observed values (). Slight improvement in the predictions at high MWN and NH2 values can be obtained with a non-linear quadratic PLS model [29]. However, the linear PLS model is very good in the target region (mid-values) and adequate in the extremes. Moreover, the control performance obtained using linear PLS model and that obtained using a non-linear quadratic PLS model, for the conditions used in this example, were found to be quite similar [29]. Estimation and model prediction assessment: One of the advantages of using PLS models for control, it is that it provides a powerful way to asses the validity of

Fig. 3. Predictions of the missing measurements made at the first decision point (35 min) using different missing data imputation methods: (-·-) expectation-maximisation (EM), (-h-) iterative-imputation (IMP), (  ) single component projection (SCP), (- -) projection to the plane using PLS (PTP) and (––) actual value.



the PLS model for trajectory estimation of the missing measurements, and for quality prediction, and it enables one to detect sensor failures, etc. because, unlike normal regression and neural network methods, it provides a model for the regressor space (x) as well as giving a prediction of the final qualities [31]. Therefore, prior to computing new control trajectories, the square prediction error (SPE) of the new vector of measurements should be computed at each decision point. This SPE provides a measure of any inconsistency between the measurements and imputed missing values for the new batch and the behaviour of the set of measurements used to develop the PLS model [15]. If the SPE is larger than a statistically determined limit [16], the quality prediction and the control computation from the PLS model should be considered to be unreliable. In this situation, it might be preferable not to recompute the MVTs at the

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553 Prediction Results

1

54

1 2 2 3 3 4

End Amine Concentration (NH2)

53

5

52

45

6

51

7 50

8 49

9

48

10

47

11

12 13 14

46

15 45 1.325

1.33

1.335

1.34

1.345

1.35

1.355

1.36

1.365

1.37

1.375 4

Number Average Molecular Weight (MWN)

x 10



Fig. 4. Observed ( ) and predicted (h) end-quality properties using PLS model.

547

control’’ region (dotted lines) was defined considering that the final product is acceptable if their predicted values lay in the specified ranges 48 6 NH2 6 50.6 and 13,463 6 MWN 6 13,590. In Fig. 5, the o’s show what happens if no control action is taken and the h’s show the end qualities obtained after control is performed. The final qualities (h’s) were obtained by rerunning the non-linear simulation model with the MVTs computed by the controller. As can be seen in this figure, the proposed control scheme corrects all batches and brings the final quality into the acceptable region. Fig. 6 shows the jacket and reactor pressure MVTs for runs 1 and 15 together with their nominal conditions. In this figure, (––) represents the MVTs computed to reject a disturbance of )10% in the initial mass of water, and (- – -) that needed to reject a disturbance of +10%. Their nominal conditions of the MVTs are indicated with (- - -).

current decision point, but simply continue to apply those from the last decision point. Control Regulatory control: At each decision time (hi ) a prediction of the final quality is made. If it is determined that control action is needed any of the control algorithms given by Eqs. (5)–(14) can be used to compute a correction, Dt, in the latent variable score space, and then the new MVTs for the remainder of the batch can be reconstructed from Eq. (21). The performance of the linear minimum variance controller algorithm (Eqs. (14) and (21) with d ¼ 1:0) is shown in Fig. 5. The final quality properties of the 15 batches shown in Fig. 4 that are affected by disturbances in the initial water concentration are shown with and without control. An ‘‘in-

Control Results 54

1

End Amine Concentration (NH2)

53

3

2

4 5

52

6

51

7 50 49

9

48

10

47

12

46

14

11

13

15 45 1.325

1.33

1.335

1.34

1.345

1.35

1.355

1.36

Number Average Molecular Weight (MWN)

1.365

1.37 4

x 10



Fig. 5. Control results: end-quality properties without control ( ); after control is taken (h) and set-point ( ).

Fig. 6. MVTs: (- - -) nominal conditions, (––) when the disturbance is )10% mass in W, and (- – -) when disturbance is +10% in W.

548

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

Note that the controller computes new MVTs that are very smooth and consistent with their behaviour during past operations. This consistency with past operation is, of course, forced to be true through use of the PLS model for MVT reconstruction (Eqs. (17) and (21)). Set-point change or new product design: In this section the performance of the control algorithm is shown in the case that a set-point change (or new product design) is desired within the region of validity of the PLS model. No disturbances in W are included in this example, but, if present, the on-line control algorithm will easily reject them as illustrated above. The desired quality properties () and those obtained by using Eqs. (14) and (21) with d ¼ 1 (h) are shown in Fig. 7 for three different setpoints. In Fig. 8 the MVTs needed to achieve such setpoints are shown. It can be seen that the performance of the algorithm in achieving the desired final quality setpoints is very good (Fig. 7), and that the MVTs computed by the controller are smooth and very consistent with the shape of the trajectories from past operation. Discussion Several practical issues may affect (to some extent) the performance of the proposed control algorithm. Some of them are briefly discussed here and more details are given in [29]. The number of latent variables is generally decided by cross-validation methods at the model building stage. It was observed that too large a number of components (with respect to that obtained by cross-validation) might promote an ill-conditioned P2 W2 inversion at the second decision point. This problem can be easily overcome by using a pseudo-inverse procedure based on singular value decomposition as detailed in [29]. For the simulation system studied no significant degradation in per-

Set-Point Change Results 1

End Amine Concentration (NH2)

54

52

50

48

46

2

44

3 42 1.29

1.3

1.31

1.32

1.33

1.34

1.35

1.36

1.37

Number Average Molecular Weight (MWN)



1.38

1.39 4

x 10

Fig. 7. Set-point change: ( ) desired, (h) achieved qualities using the control algorithm (Eq. (14) and (21)) and ( ) nominal operating point.

Fig. 8. MVTs for set point change: the number indicates the set-point change shown in Fig. 7, and (- - -) the nominal MVT.

formance is obtained by using a different number of PLS components. The influence of using different missing data imputation algorithms was also studied. All the algorithms give adequate control performance. Those based on EM, IMP and PTP perform slightly better than the one in which SCP was used. In the previous examples, a total of 45 batches (30 with a movement in the MVTs at the two decision points) were used for model identification. However, adequate control performance (all test batches falling inside the ‘‘in-control’’ region of Fig. 5) was achieved using as few as 15 batches (10 in which some experiment in the MVTs was performed). This illustrates that the data requirements for PLS model building are modest. However, if the model has been identified using very limited or uninformative batch data-sets (as those arising from only historical data), batch-to-batch model parameter updating can be performed at the end of each new completed batch to improve the quality of the

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

model parameter estimates, prediction and control for the upcoming batch [12,29]. To assess the impact of measurement noise on the performance of the algorithm, different levels of random noise were added to the on-line measurements of reactor temperature (Tr ) and venting (v). It was found that adequate control performance (test batches falling inside the ‘‘in-control’’ region of Fig. 5) was achieved with noise levels up to 35% in the temperature and venting rate. The noise level here represents the percentage of the noise variance over the true variations of the temperature and venting rate changes observed in the training set. The 35% noise level approximately represents one standard deviation in temperature of 2 K and venting rate of 42 g/s (see [29] for details). For larger levels of noise (50%, for example) the control performance is degraded to some extent because the random error added to the measurements becomes quite large when compared with the true variations in the MVs. A no-control region that reflects the impact of these measurement noises may be obtained by propagating such measurement errors with the PLS model as suggested in [9]. This would prevent control actions from being implemented based solely on the uncertainty arising from noise. Finally, the control methodology outlined in Section 2 can be easily extended to cases in which a non-linear PLS model and control is needed. This is achieved by simply modifying Eq. (12) (case 3, dimðDtÞ > dimðysp Þ) to take into account the non-linear nature of the PLS algorithm. For example, in the case of a quadratic PLS model, Eq. (12) can be restated as: min |{z}

DtT Dt

st

yTsp ¼ uT QT

ð22Þ

Dtðhi Þ

where uT ¼ b1 þ b2 ðDt þ ^tpresent ÞT þ b3 ðDtT þ ^tTpresent Þ2 .

549

This equation can be easily solved for Dt using quadratic programming (or non-linear least squares in the case of no constraints), and the MVTs can be reconstructed in the same way as described in Section 2.4. From the simulation study, the control performance of the quadratic PLS model is quite similar to that obtained using the linear PLS model [29] (and therefore results are not shown). This is not surprising because, in the region under study, the process is only slightly non-linear. However, if larger disturbances affect the process a non-linear PLS approach may be better suited.

3.2. Case study 2. Feasibility study on industrial data for an emulsion polymerisation process Data: In this feasibility case study, industrial data for an emulsion polymerisation processes is used. The original data set consists of 53 batches obtained from an experimental design in which the initial conditions and/ or process variable trajectories were altered. No intermediate quality measurements were available during the reaction. However, final product physical properties (FP) and final product quality properties (FQ) are available at the end-of the process for most of the batches. Fig. 9(a) shows the actual process variable trajectories that comprise the training data set (X), while Fig. 9(b) shows the 6 quality properties (Y matrix), corresponding to these batches. In Fig. 9(a), it can be noticed that (i) since the batches in the process were of unequal duration, alignment of the trajectories was accomplished using the reaction extent as an indicator variable [15,21] (every interval represents a 0.5% increase in the reaction extent), and that (ii) some of the trajectories contain a noticeable level of noise. It was decided not to perform

Fig. 9. (a) Original process variable trajectories (every interval represents 0.5% of reaction extent), and (b) original quality properties.

550

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

any pre-treatment such as filtering or smoothing on the process trajectories in order to test the performance of the prediction and control algorithm under this situation. It can also be seen in Fig. 9(a) that FP-1 and FP-2 are highly correlated therefore, to avoid an ill condition matrix inversion in the control computation stage, FP-2 was removed and only five end quality properties controlled. Removing FP-2 poses no problem since by controlling FP-1 and the other quality variables we are controlling FP-2 indirectly. Alternatively, we can perform PCA on the quality property matrix (Y) and control the corresponding principal components instead of the actual properties. For property reasons no further details will be given here regarding the nature of the process trajectories, initial conditions or product specifications. From the original data, 49 batches were used as a training data set, while four batches were used as testing set. These four batches were selected to span different regions of the space far from the origin as can be seen in Fig. 10. In this figure the projection of all batches in the first two PLS dimensions (t1 –t2 ) is shown. Batches 6, 12, 16 and 46 were removed from the dataset and used as test data. The 8 process variable trajectories are manipulated variables and each one of them is segmented in 200 intervals (every interval represents 0.5% of reaction extent). Therefore the data matrix used for model building consists of segmented MVTs [X] and initial conditions [Z] (regressor matrices), and the matrix of five physical and quality properties [Y]. The identified PLS model consists of five latent variables (obtained by cross-validation) that fits 76.8% of the X space and 69.9% of the Y space. Based on cross-validation, 51.7% of the Y space can be predicted. Control objectives: The batch data in this study was the result of open-loop batch runs collected under different initial conditions and different MVTs. There was no possibility of implementing the resulting controller on the batches. Therefore, this data is simply used to test the feasibility of the prediction and control algorithms.

60

18

40 5 20

7

t [2]

9 8

0

12

3 11 1

4 2 14 13

6 2115

16

17

10 -20

19 36 29 23 38 26 37 355120 22 39 24 25 40 28 27 3253 30 34 44 31 524142 5043 47334948 45 46

One of the existing batch runs is taken as the nominal conditions and the final physical and quality variables (y) measured from it selected as the targets (set-points). Others batch runs with different initial conditions and different MVTs are then selected as initial disturbance conditions for a new batch. If no corrective action is taken to adjust the MVTs then the batch will follow the actual MVTs implemented throughout its duration, and the final quality (y) will be the measured values for that batch. Control is to be applied after a batch has reached 10% of completion (based on reaction extent). Direct evaluation of the controller is not possible, but indirect validation can be obtained by comparing how close the recomputed MVTs follow the nominal MVTs from 10% of reaction extent until the end of the batch. Since the first 10% of the history of the new batch is different from the nominal MVTs, then to achieve the desired final qualities (qualities of the nominal batch), one should not expect the recomputed MVTs to exactly follow those for the nominal batch, but they should be close to them. Notice that if the control algorithm is actually implemented, it would pose no problem to recompute the MVTs at several decision points and not only at one as shown here. Prediction: To evaluate the performance of the PLS missing data algorithms, the total percent relative RMSE for all the qualities properties (5 in this study) is shown in Table 1 over the k ¼ 4 batches that compose the testing data set: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 !1 u k 5 uX yij
^yij 2 X 1 1 t A  100 %RMSE ¼ @ 5 i¼1 k j¼1 yij where yij is the i observed end-quality property for the j batch and ^yij its predicted value.As an illustration of the missing measurement reconstruction (at 10% of reaction extent using the EM approach), Fig. 11 is shown for batch 12, where it can be noticed that the trajectory estimation is satisfactory in spite of the high level of noise. Control: As an illustration of the control performance using the proposed scheme (Eqs. (10) and (21) with d ¼ 1.0), results for one testing batch (batch 12) are shown. Fig. 12 shows the measured final values of the y variables () for the batch when no control was taken, their predicted values at 10% of completion if no control were taken ( ), the target values (h), and the expected quality properties obtained if control action were performed ( ). Since a minimum variance strategy was used

-40 -60

-50

-40

-30

-20

-10

0 t [1]

10

20

30

40

50

60

Fig. 10. t1
t2 PLS space for the batches used in the training data set. Batches 6, 12, 16, and 46 were removed from the original data set and used as test data.

Table 1 Performance of missing data algorithms for prediction: total percent relative RMSE for all five end quality properties Algorithm

EM

IMP

SCP

PTP–PLS

%RMSE

8.0

7.3

9.8

6.8

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

551

Fig. 11. Performance of the missing data algorithm for reconstruction of process measurements. The prediction is performed at 10% of reaction extent (every interval represents 0.5% of reaction extent): (  ) estimated trajectory using the EM algorithm and (––) observed trajectories (scaled units).

FQ-1

Quality Property Value

9

FQ-2

1.4

8

1.2

10

7

1

9

6

0.8

8

5

0.6

7

4

0

1

2

0.4

0

2

6

4

0.7

3.5

0.65

3

0.6

2.5

0.55

0

1

0

1 5

FP-1

0.75

1

FQ-3

11

2

2

x 10

0

2

FP-3

1

2

Fig. 12. Control results (control action taken at 10% of completion of the batch). Target (h), predicted qualities ( ), observed values if no control action is taken ( ) and expected quality properties if control action were performed ( ).



(Eq. (10) and (21)), the values of the expected end quality properties resulting from the control algorithm will match their targets ( ), (since these values were computed using simply the PLS model with the imputed

MVT adjustments obtained from model inversion using the same PLS model.) A better way to evaluate the reasonableness of the control is to inspect the MVTs obtained from the control algorithm. Fig. 13 shows

552

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

Fig. 13. MVTs (computed at 10% of reaction extent from the beginning of the process): (  ) nominal conditions; (- - -) current trajectories that would give ‘‘out-of-control’’ qualities and (––) MVTs obtained from the control algorithm (equation (10) and (21), with d ¼ 1:0).

nominal trajectories (  ), the current trajectories that would give ‘‘out-of-control’’ qualities (- - -) and the MVTs obtained from the control algorithm (––) (at 10% of reaction extent) that would drive the predicted physical and quality properties to the desired targets. In this figure, notice that MVTs obtained from the control algorithm after 10% of completion are quite close to their nominal conditions and exhibit the desired shapes. It seems reasonable to assume that if these new trajectories were to be implemented, they would drive the process closer to the desired end-quality values, simply because the new MVTs are much closer to the nominal conditions than those when no control is performed. Note that they should not match the nominal trajectories exactly because they must also compensate for the first 10% of the batch being run at the wrong conditions. Furthermore, since the trajectories are highly correlated with one another, there are various trade-off among the MVTs that might give quite similar final quality values. In summary, although the control could not actually be tested, these results indicate that the controller is behaving very much as one might expect and are providing the incentive for its implementation.

that recomputes, on-line, the entire remaining trajectories for the MVs at several decision points. In spite of the fact that the resulting controller solves for the high dimensional MVTs, the control algorithm involves solving for only a small number of latent variables in the reduced dimensional space of a PLS model. The high dimensional MVTs are then solved by inverting the PLS model. The only requirement of this approach (as with any other control algorithm that recomputes the MVTs) is that the lower level control scheme can accept and track the computed modified trajectories. The strategy uses empirical PLS models identified from historical data and a few complementary experiments. The algorithm is illustrated using a simulated condensation polymerisation process and data obtained from an industrial emulsion polymerisation setting. Since smooth and continuous MVTs can be obtained, the approach seems well suited for use in processes and mechanical systems (robotics) where such smooth changes in the MVs are desirable. The methodology would also be well suited to the control of transitions of continuous processes.

Acknowledgements 4. Conclusions A novel control strategy for final product quality control in batch and semi-batch processes is proposed

J. Flores-Cerrillo thanks McMaster University and SEP for financial support and to Dr. Russell S. A. for kindly providing us with his condensation polymerisation simulator.

J. Flores-Cerrillo, J.F. MacGregor / Journal of Process Control 14 (2004) 539–553

References [1] D.J. Kozub, J.F. MacGregor, Feedback control of polymer quality on semi-batch copolymerization reactors, Chem. Eng. Sci. 47 (4) (1992) 929–942. [2] C. Kravaris, M. Soroush, Synthesis of multivariate nonlinear controllers by input/output linearization, AIChE J. 36 (1990) 249– 264. [3] C. Kravaris, R.A. Wright, J.F. Carrier, Nonlinear controllers for trajectory tracking in batch processes, Comp. Chem. Eng. 13 (1989) 73–82. [4] T.J. Crowley, K.Y. Choi, Experimental studies on optimal molecular weight control in a batch-free radical polymerisation process, Chem. Eng. Sci. 53 (15) (1998) 2769–2790. [5] D. Ruppen, D. Bonvin, D.W.T. Rippin, Implementation of adaptive optimal operation for a semi-batch reaction system, Comp. Chem. Eng. 22 (1/2) (1997) 185–199. [6] D. Bonvin, B. Srinivasan, D. Ruppen, Dynamic optimization in the batch chemical industry, in: Proceedings of the CPC-VI Conference, AIChE Symp. Ser. 326 (98) (2002) 255–273. [7] B. Srinivasan, S. Palanki, D. Bonvin, Dynamic optimization of batch processes: I. Characterization of the nominal solution, Comp. Chem. Eng. 27 (1) (2003) 1–26. [8] B. Srinivasan, D. Bonvin, E. Visser, S. Palanki, Dynamic optimization of batch processes: II. Role of measurements in handling uncertainty, Comp. Chem. Eng. 27 (1) (2003) 27–44. [9] Y. Yabuki, J.F. MacGregor, Product quality control in semibatch reactors using midcourse correction policies, Ind. Eng. Chem. Res. 36 (1997) 1268–1275. [10] Y. Yabuki, T. Nagasawa, J.F. MacGregor, An industrial experience with product quality control in semi-batch processes, Comp. Chem. Eng. 24 (2000) 585–590. [11] J. Flores-Cerrillo, J.F. MacGregor, Control of particle size distributions in emulsion semi-batch polymerization using midcourse correction policies, Ind. Eng. Chem. Res. 41 (2002) 1805– 1814. [12] J. Flores-Cerrillo, J.F. MacGregor, Within-batch and batch-tobatch inferential-adaptive control of batch reactors: A PLS approach, Ind. Eng. Chem. Res. 42 (2003) 3334–3345. [13] S.I. Chin, S.K. Lee, J.H. Lee, A technique for integrated quality control, profile control, and constraint handling for batch processes, Ind. Eng. Chem. Res. 39 (2000) 693–705. [14] S.A. Russell, P. Kesavan, J.H. Lee, B.A. Ogunnaike, Recursive data-based prediction and control of batch product quality, AIChE J. 44 (1998) 2442–2458. [15] P. Nomikos, J.F. MacGregor, Multivariate SPC charts for monitoring batch processes, Technometrics 37 (1) (1995) 41–59.

553

[16] P. Nomikos, J.F. MacGregor, Monitoring batch processes using multiway principal component analysis, AIChE J. 40 (8) (1994) 1361–1375. [17] P. Nomikos, J.F. MacGregor, Multi-way partial least squares in monitoring batch processes, Chemometr. Intell. Lab. Syst. 30 (1995) 97–108. [18] G. Chen, T.J. McAvoy, Process control utilizing data based multivariate statistical models, Can. J. Chem. Eng. 74 (1996) 1010–1024. [19] G. Chen, T.J. McAvoy, M.J. Piovoso, A multivariate statistical controller for on-line quality improvement, J. Proc. Cont. 8 (2) (1998) 139–149. [20] T.J. McAvoy, Model predictive statistical process control of chemical plants, in: Proceedings of the 15th IFAC World Congress, Barcelona, Spain, 2002. [21] A. Kassidas, J.F. MacGregor, P.A. Taylor, Synchronization of batch trajectories using dynamic time warping, AIChE J. 44 (1998) 864–875. [22] T. Kourti, P. Nomikos, J.F. MacGregor, Analysis, monitoring and fault diagnosis of batch processes using multi-block and multiway PLS, J. Process Control 5 (1995) 277–284. [23] T. Kourti, J. Lee, J.F. MacGregor, Experiences with industrial applications of projection methods for multivariate statistical process control, Comp. Chem. Eng. 20 (1996) S745–S750. [24] D. Neogi, C.E. Schlags, Multivariate statistical analysis of an emulsion batch process, Ind. Eng. Chem. Res. 37 (1998) 3971–3979. [25] P. Nelson, J.F. MacGregor, P.A. Taylor, Missing data methods in PCA and PLS: score calculations with incomplete observations, Chemometr. Intell. Lab. Syst. 35 (1996) 45–65. [26] F. Arteaga, A. Ferrer, Dealing with missing data in MSPC: several methods, different interpretations, some examples, J. Chemometr. 16 (2002) 408–418. [27] S. Garc
ıa Mu~ noz, J.F. MacGregor, T. Kourti, Model predictive monitoring of batch processes, Ind. Eng. Chem. Res., submitted for publication. [28] J.M. Jaeckle, J.F. MacGregor, Product design through multivariate statistical analysis of process data, AIChE J. 44 (1998) 1105– 1118. [29] J. Flores-Cerrillo, Quality control for batch processes using multivariate latent variable methods, Ph.D. thesis, McMaster University, Hamilton, ON, Canada, 2003. [30] S.A. Russell, D.G. Robertson, J.H. Lee, B.A. Ogunnaike, Control of product quality for batch nylon 6, 6 autoclaves, Chem. Eng. Sci. 53 (1998) 3685–3702. [31] A. Burnham, R. Viveros, J.F. MacGregor, Latent variable multivariate regression modeling, Chemometr. Intell. Lab. Syst. 48 (2) (1999) 167–180.