4202
Ind. Eng. Chem. Res. 2008, 47, 4202–4208
Optimization of Batch Operating Policies. Part II. Incorporating Process Constraints and Industrial Applications Salvador Garcı´a-Mun˜oz,†,‡ John F. MacGregor,*,† Debashis Neogi,§ Bruce E. Latshaw,§ and Sanjay Mehta§ Department of Chemical Engineering, McMaster UniVersity, 1280 Main Street West, Hamilton, ON, L8S 4L7 Canada, and Air Products and Chemicals, 7201 Hamilton BouleVard, Allentown, PennsylVania 18195
In the first part of this series [Ind. Eng. Chem. Res. 2006, 45, 7856-7866], data-driven approaches, based on partial least squares (PLS) models built from historical batch data, were used to find optimal batch operating trajectories that would yield a desired vector of final product quality attributes. The method allowed for the inclusion of univariate and multivariate constraints on the set of desired final product quality attributes and presented approaches for handling multiple solutions. In this paper, the technology is further extended to include constraints in the process operating trajectories themselves. The methodology is successfully applied to an industrial batch polymerization process where the batch trajectories are designed to achieve specific properties of the final polymer while consuming the minimal amount of time for the batch run. 1. Introduction Batch processes play a key role in the manufacturing of high value added chemicals such as pharmaceuticals and polymer resins. Process systems engineering (PSE) has played a crucial role in developing techniques to improve the operation for batch units. At the design phase where the geometry and nominal operation is determined, fundamental models are traditionally used to aid the design exercise. Once a class of products is in production, any improvements or optimization of the product quality or process operation can use empirical models built from the accumulated database. Multivariate latent variable models (LVM) have proven to be efficient at modeling batch process operation2 and the sources of uncertainty that enter the system. 1.1. Background. The redesign of process operating conditions to achieve a new set of final product quality attributes, based on partial least squares (PLS) models built from historical data, was first discussed by Jaeckle and MacGregor.3,4 In those works the objective was to use historical databases on process operation over a range of product grades to find new operating conditions that might achieve a desired set of new product properties. This was referred to as a product design problem. The solution was provided by an inversion of the linear or nonlinear PLS models. The first paper in this series1 extended the approach to the estimation of complete batch trajectories and considered optimization based approaches to handle multiple solutions and constraints in the final product properties. An alternative approach, proposed to introduce more information about the derivatives of the trajectories into the multiway model, was shown to provide smoother solutions and eliminate the multiplicity of solutions. 1.2. Contributions of this Work. The way that the variable trajectories were estimated in the first paper1 does not allow for the introduction of constraints in the process variable trajectories into the optimal design criterion. The focus of this second paper is to reformulate the optimization to include * To whom correspondce should be addressed. Tel.: (905) 525-9140 ext 24951. Fax: (905) 521-1350. E-mail:
[email protected]. † McMaster University. ‡ Current address: Pfizer Global Research and Development, 445 Eastern Point Road, Groton, CT 06340. E-mail:
[email protected]. § Air Products and Chemicals.
additional objectives and constraints on the batch operation into the design, and to present results on an industrial application of the method. The industrial problem involves the design of new operating trajectories and initial conditions (recipes) for a batch polymerization process to achieve a desired set of end product properties while using the shortest possible batch duration. The development of several new polymer grades is considered and some of the solutions are implemented in an industrial scale pilot plant. 2. PLS Modeling of Batch Processes This paper follows the same notation as in past work1–7 on multivariate modeling of batch processes. A data set from a batch process consists of I batches, and measurements on J variables taken throughout the duration of the batch (variable trajectories). Due to variations in the time duration of the batches, there may be different numbers of observations throughout each batch. However, after alignment of the batch data using a maturity or evolution index variable6 (or other approach5) one is left with a fixed number of observations at K values of the maturity variable throughout each batch. After unfolding the three-way batch data array X as proposed by Nomikos and MacGregor,6,7 the two-way matrix of variable
Figure 1. Structure of unfolded and aligned batch data.
Figure 2. Data structure with PLS score matrix T and new vectors corresponding to the desired trajectory optimization solution.
10.1021/ie071437j CCC: $40.75 2008 American Chemical Society Published on Web 05/15/2008
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4203
Figure 3. Raw process trajectories for the semibatch process.
trajectories (X) consists of I rows and JK columns. As discussed in part I,1 there are often advantages to augmenting the X matrix to include the smoothed derivatives of all or some of the process variable trajectories (DX). Furthermore, at the beginning of the run, there might be a set of L initial conditions (raw material properties, ingredient charges, or prior processing history of the materials, etc.) that are recorded for each run. At the end of the run, the quality and yield of the final product is characterized by a set of M properties. The initial conditions will form an I × L matrix referred as to Z, and the set of final product quality characteristics form an I × M matrix referred to as Y. The data structure is shown in Figure 1. The combined data blocks Z, X, and DX (with suitable block scaling) will be denoted as the augmented X-block, XA. Each row represents the data available on one batch. A PLS model built using data from the I batches will have the following eqs where T is a matrix with a columns of scores for the a statistically significant latent variables (Figure 2): XA ) TPT + E
(1)
Y ) TQT + F
(2)
T ) XAW*
(3)
T T ydes ˆ ) τnew Q
(4)
The latent variable score matrix T is determined from the PLS algorithm as a linear combination (eq 3) of the initial conditions, the unfolded trajectories, and trajectory derivatives as defined by the PLS weighting matrix W*.8 P and Q are loading matrices obtained from the PLS algorithm that show how XA and Y, respectively, are related to the scores T.8 In the T batch trajectory optimization problem considered here, yˆdes is the projection onto the PLS model of the vector of desired final T batch quality attributes ydes for which a new set of batch T T trajectories and initial conditions are needed (xanew ). τnew is the
vector of PLS scores corresponding to the new desired quality T attributes yˆdes . These vectors represent a new row of process conditions, product qualities and scores corresponding to the desired new batch as illustrated in Figure 2. 3. Trajectory Estimation In earlier papers,1,3 the approach focused in the estimation T of the PLS score vector (τnew ) that corresponded to the desired set of characteristics in the final product. Once the score vector was determined, the trajectories (and the initial conditions if present) were computed by simple projection onto the plane describing the operating space: xanew ) Pτnew
(5)
T The part of xanew corresponding to the variable trajectories and the initial conditions provided the solution. In the work of T Jaeckle and MacGregor,3 the solution to the PLS scores τnew T that correspond to the new desired product quality yˆdes was obtained by the pseudo inverse of the PLS model eq (4):
τnew ) (QTQ)-1QTydes
(6)
This solution allowed for no constraints on the elements of τTnew (such as how far the new solution could be from the existing data) or for one sided constraints, etc., on individual elements T of yˆdes . In the work of Garcia et al.,1 the inversion of the PLS model to obtain the PLS scores τTnew was posed as an optimization problem, (since this estimate will be used in a second T T optimization, the new notation will refer to τnew as τˆ xnew ):
{
(∑ )} A
min (ydes - Qτˆxnew)TG1(ydes - Qτˆxnew) + F ˆτxnew
ˆτxnew,a2
sa2 s.t. BQτˆxnew < b a)1
(7)
4204 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008
Figure 4. Aligned process trajectories and time usage for the semibatch polymerization process.
Where G1 is a square matrix with the weights given to the multiple elements of ydes in the main diagonal; sa2 is the eigenvalue of principal component a; and F is a weight given to the Mahalanobis distance of the solution. In this way, the T distance of the solution τˆ xnew from the center of the existing data is penalized by the second term in the quadratic cost function (a penalty on the Mahalanobis distance of the scores T from the origin), and inequality constraints on elements of yˆdes were addressed using a quadratic programming (QP) solution. T In the current paper, the PLS scores τˆ xnew that correspond to T the desired quality vector ydes are again determined by the optimization (7). However, constraints on the process variable trajectories are also considered and are either imposed directly T as constraints on the elements of xanew or as added terms T T involving xanew in a second optimization for xanew (rather than using the unconstrained solution (5)). Furthermore, in order to T satisfy these constraints on xanew some small extrapolation of the model (nonzero residual squared prediction error (SPE)) will be allowed. The general idea is to tolerate a certain amount of model mismatch (nonzero SPE) between the estimate of the design and its prediction from the PLS model and moderate extrapolations in the score space to allow for achieving constraints on xaTnew. These new exploratory designs may contain slight distortions in the trajectories in order to satisfy the optimization, but as long as the SPE is kept within a conservative limit and if the derivative augmented multiway PLS model approach is used to force smoothness, the solutions should be quite consistent with past operating conditions. The following objective function for the second stage optimization can then be used (as an alternative to eq 5) to obtain T a solution for the new trajectories and initial conditions xanew : min {(τ T xanew
T τxnew)TG2(τxnew - ˆτxnew) + (xanew xnew - ˆ T T T T )TΛ(xanew -ˆ xanew ) + ηxanew } xˆanew
(8)
T In this objective function, xanew is the solution vector of initial conditions and process and derivative trajectories for the new T batch product; xˆanew is the projection of the solution onto the PLS model (the two solutions will differ slightly if the SPE is nonzero); τˆ xnew is the PLS score vector solution corresponding T to ydes coming from the first optimization in eq 7, and τxnew ) T T W* xanew is the score vector corresponding to the solution xaTnew in eq 8. The first term in (8) tries to minimize the error between τxnew ) W*Txanew and the desired τˆ xnew from (7). The second T term in (8) is the SPE or distance of the solution xanew from its projection onto the PLS model plane. Minimizing this will allow only small extrapolations in the solution for xaTnew from the PLS model. The last term in (8) penalizes large values of any T trajectory in xanew . A particular concern in the industrial application of the methodology in this paper is obtaining a T solution that not only yields a yˆdes in the proper region but also minimizes the batch duration needed to achieve this. It is shown later that, upon alignment of the trajectories, a new variable trajectory, “time usage”,5 is added to the vector of process variable trajectories. By specifying nonzero elements in the weighting vector η only for the elements of xaTnew corresponding to time usage, then this will allow for minimization of batch time at every stage throughout the batch. Other functions of the process variables like a minimal consumption of a raw material can be introduced within the same framework; the only factor to adjust are the values given to the elements of the η weighting vector for the chosen variable trajectory in eq 8. If, in addition, certain process variable trajectories had to follow prespecified target trajectories as closely as possible, then an additional term (xanew - xatarget)TΩ(xanew - xatarget) could be added. G2, Λ, and η in eq 8 are weighting matrices and a weighting vector that determine the relative importance of the various terms in the objective function. Some suggestions for these are discussed in part I.1
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4205
Figure 5. Several designs for minimal time: (-) solution 1; (- - -) solution 2; ( · · · ) solution 3; (- · -) solution 4; (s) solution 5.
The objective function (8) can be expressed in a standard QP form.
{
1 T Hoxanew + foTxanew min xanew xanew 2
}
where Ho ) W*G2W*T + Λ - ΛPW*T - W*PTΛ + W*PTΛPW*T T foT ) η - τxnew G2W*T
(9)
The loadings (W* and P) to build Ho and fo should be those from the PLS model (Figure 2) where the initial conditions (Z) and the variable trajectory blocks (X and DX) are considered in the XA space (e.g., by scaling each block to a chosen variance and using this block scaled [Z X DX] as XA in the PLS model). A data set with batches having varying time durations is required in order to perform any optimization in the time dimension. Without varying batch times the data set will contain no information on the impact that varying batch time will on the product quality Y. These variable duration batches then must first be aligned5,9 before the PLS model is fitted. Once the alignment operation is done, a time-usage variable trajectory5 will be available. The time-usage trajectory for each batch
Table 1. Square Prediction Error for the Five Different Optimal Solutions and the 95% Limit from the Training Data Set solution
SPE
SPE 95% CI
1 2 3 4 5
4.24E-019 0.01 0.04 0.15 0.6
6.88 6.88 6.88 6.88 6.88
contains the actual clock time taken to reach each value of the maturity variable in the aligned batches.5 This and its derivative are then added to the unfolded batch matrix XA as additional trajectory variables. The new time-usage trajectory is the key trajectory to be minimized in order to achieve a design with minimal time duration. This minimization is accomplished in eq 8 by assigning positive values to the η elements corresponding to all the elements in the time-usage trajectory. The solution to this final formulation should provide a design that will yield the desired quality, with a low usage of time, while keeping a conservative Mahalanobis distance from the existing training data, and a conservative orthogonal distance (SPE) to the model plane. This approach is illustrated in the following section with the use of an industrial data set from a polymerization process.
4206 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 Table 2. New Desired Properties property
new grade 1
new grade 2
Y1 Y2 Y3 Y4 Y5 Y6 Y7
0.67 0.75 62.2 2.35 3000 15.47 0.87
0.77 0.991 64 6.5 2800 15.38 0.33
4. Application to an Industrial Polymerization Process To evaluate the potential of this methodology, studies were carried out in an industrial pilot plant semibatch polymerization reactor at Air Products and Chemicals. The data set available consisted of a number of batches run with a variety of initial conditions and different operating trajectories. Some of these were collected from a designed experiment. The initial condition matrix (Z) consisted of 11 variables, including the initial charges of various monomers, stabilizers, initiators, etc. These charges were ratioed to the total charge of a base monomer to ensure consistency among batches. Eight process variables were measured throughout the course of the reaction (process variable trajectories, X) as shown in Figure 3. These consisted of a measure of conversion, monomer and initiator addition profiles, and reactor physical conditions, denoted as X1-X8. In addition, seven final polymer properties and six end-use quality variables (Y) were measured on the product produced from each batch. The polymer properties consisted of typical polymer characteristics such as molecular weight and glass transition temperature (Tg). The end use quality variables consisted of market oriented applications tests used to characterize the polymer’s performance. Note that the duration of all the batch runs varies considerably. This will necessitate that the batches be aligned for PLS analysis, but this time variation will be preserved in the model and will also provide the information needed to find optimal trajectories with minimum batch duration. Prior to the alignment, a low-pass median filter is applied to smooth out some noisy peaks in the data. Using the maturity (or indicator) variable approach,5 the trajectories are aligned with respect to the conversion. Sampling every 0.1% for the first percent, and every 0.5% for the rest of the batch, for a total of 208 samples per batch (this multirate sampling scheme was used simply to better capture the start of the reaction and adds no complexity to the PLS modeling). The aligned trajectories are shown in Figure 4. A time usage trajectory, computed as in Garcia et al.,5 has been added to the aligned data set (Figure 4) to replace conversion since conversion is now the maturity variable (abscissa of the plots). The time usage trajectories show the actual clock time needed to reach any given level of conversion during each batch. It is shows the speed at which each batch is proceeding with respect to other batches at all conversion levels throughout the batch duration. This time usage variable, constructed from the alignment, preserves the critical batch timing information in the subsequent analysis (and allows one to reconstruct the original nonaligned trajectories). We now consider two studies using the models built from these data. In the first we present a complete application of the methodology to the design of the optimal operating policies for a new batch aimed at achieving a product having the 13 final polymer attributes within a specified region and having a minimum batch time. In the second study we report on the results obtained from applying the method in the pilot plant reactor for the development of two new products.
4.1. Illustration of the Methodology Applied to the Development of a New Product. The objective in this section is to design the batch operating policies for a new product that satisfies the following requirements on the 13 product attributes: 1. y1 to y6 and y8 are to be kept within the normal range. 2. y7 ) y7des and y9 ) y9des. 3. y10 < y10const, y11 < y11const, y12 < y12const, and y13 < y13const. In addition, the batch duration was to be minimized, subject to satisfying the product attribute conditions above. A PLS model was built using the aligned trajectory data (shown in Figure 4) and using the augmented XA matrix shown in Figure 2 in which both the initial condition and trajectory derivatives are included. The equality and inequality constraints shown above (requirements 1-3) are included in the first optimization objective function (eq 7) that is solved to estimate τˆ xnew. Once the τˆ xnew is estimated, the trajectory design problem is solved by minimizing the second optimization in eq 8 for xaTnew using its quadratic programming form (eq 9). In this study, no constraint was placed on the Mahalanobis distance in eq 7 (i.e., the weight F ) 0) since derivatives were included in the augmented XA matrix, thereby eliminating the presence of multiple solutions.1 Several optimal trajectories solutions are shown in Figure 5. Besides the usage of time, the solutions differ in their orthogonal distance to the plane (Table 1), which is still being kept at a very conservative value through the penalty weight matrix Λ in eq 8. In solution 1 (Figure 5), no constraints are set on the batch time usage or any other process variable. Its solution can be arrived at through solving eq 7 for τˆ xnew and, then, obtaining the process trajectories xaTnew either though direct inversion of the PLS model of the XA space (eq 5), or through solution of the second optimization in eq 8 with no constraints on the process trajectories (η ) 0). Solutions 2-5 show the results of applying a successively increasing constraint on the time usage trajectory for the new batch. They are obtained by the two-step optimization of solving the optimization in eq 7 followed by the optimization in eq 8 with increasing weight on time usage trajectory (i.e., increasing the magnitudes of the elements of η corresponding to the time usage variable trajectorysweights on all the other process variable trajectories are set to zero.). As the time usage is penalized more heavily, the orthogonal distance T from the model plane of the optimal process trajectories xanew increases slowly (Table 1). This is expected, since the solution is beginning to deviate slightly from the model in order to achieve shorter batch times. However, the SPE deviations are low (well below the 95% confidence limit for the training data set), and the solutions still obey the dominant covariance structure of the variables and derivatives as captured in the PLS model. This is evident from the similarities of the successive solutions shown in Figure 5. The fact that all the trajectories are being adjusted simultaneously shows that the solution is not just focusing on lowering the time usage trajectory alone, but performing adjustments on the entire set of process variable and derivative trajectories such that solution still satisfies the requirements for the product quality attributes as predicted by the PLS model. For comparison purposes, all of the five solutions were forced to project to the same desired point on the latent variable space (τxnew ) τˆ xnew) by adjusting G2 along with η in eq 8 (the G2 weighting matrix has to be incremented along with η to keep τxnew ) τˆ xnew). Although this is not generally necessary, it was done here to illustrate the effect on the residuals (SPE) when the solution for the final product is kept constant at the desired
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4207
Figure 6. Estimated and executed trajectories for five manipulated and three measured variables for grade 1.
Figure 7. Estimated and executed trajectories for five manipulated and three measured variables for grade 2. Table 3. Desired and Obtained Properties of the Two New Polymers grade 1 desired grade 1 obtained grade 2 desired grade 2 obtained Y1 Y2 Y3 Y4 Y5 Y6 Y7
0.68 0.75 62.20 2.36 3000.00 15.48 0.87
0.69 0.76 57.50 2.35 3120.00 16.09 0.89
0.77 0.991 64 6.5 2800 15.38 0.33
0.77 0.991 61.8 5.08 2875 15.38 0.33
quality (ydes ) Qτˆ xnew) for all solutions. The residuals (SPE) for all five solutions are summarized in Table 1. In this illustrative industrial example, the time usage5 information from the alignment is a key concept in order to find minimum batch time designs, but other functions of the process variables like a minimal consumption of a raw material can be introduced within the same framework; the only factor to adjust are the values given to the η weighting vector in eq 8.
The solution (if close enough to the plane) is guaranteed to have the same correlation structure as the historical data and therefore should not show any abnormal behavior that may make the design infeasible from a practical point of view. 4.2. Pilot Plant Results for the Design of Two New Products. The second study involved the design two new products within the region of the products produced by the batches in Figure 3. The PLS model with an augmented XA matrix including trajectory derivatives and the batch initial conditions was again used. The products had to satisfy the requirements listed in Table 2, while using a given set of initial conditions and using the least amount of time to complete the batch run. To first check if there is evidence that these product properties can be achieved within the current batch process and within its past operating region, the two sets of y values in Table 2 are projected onto the plane of a PCA model built from all the Y
4208 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008
data obtained from the batch training runs (historical plus DOE) shown in Figure 3. If the new grades of polymer fall close to this plane, there is strong evidence that they are feasible. For these cases the squared prediction error (SPE) for the first new grade is within the 99% confidence interval (CI) for the SPE from the PCA model of the past y’s and the SPE for the second new grade is slightly above the 99% limit. The estimation must therefore be done with caution since the desired new grades are not entirely coplanar with the previous ones. The equality constraints for the final properties are then included into the objective function (7) to estimate τˆ xnew. Once τˆ xnew is estimated, the trajectory design problem is solved by minimizing eq 8. Once a set of trajectories is estimated for each grade of new polymer, it is up to the manufacturing facility to implement them. To make this easier for the operators, it was suggested that the decomposition of the optimal trajectories be done by taking into account the steps and events defined in the batch recipe. Doing so will enable the optimization to provide results that are straightforward to implement. For this example, the pilot plant operators programmed a recipe that would resemble the estimated trajectories. In this particular case, the recipes are implemented thru a PLC system. Figures 6 and 7 illustrate the estimated and implemented trajectories. Table 3 lists the properties of the new desired and obtained grades of polymer, as reported by the analytical test on the final product yielded by the two executed batches. The obtained properties are considered to be very close to those desired in both cases. According to the operators, the new batch runs are also between 20 and 30 min faster than expected from their experience. A formal process development exercise by traditional means and its implementation would be required to prove this last point. 5. Conclusions A method for the optimization of batch operating policies is proposed based on multivariate partial least squares (PLS) models built using data collected from historical batch runs. The objective of the optimization is to find the batch initial conditions and manipulated variable trajectories that can achieve a final product having specified properties and that can be
achieved in a minimum batch time. The optimization framework proposed allows for the inclusion of constraints on the batch trajectories as well as on the final product attributes. The resulting solutions are guaranteed to be consistent with past operating policies and constraints since they are forced to be consistent with the PLS model built from past data. This is both a strength and weakness of the approach. It ensures that one obtains feasible solutions that can easily be implemented within the current control and operation scheme. However, it cannot design operating policies for products that are very different from the past, nor use operating trajectories and recipes that are unlike anything used in the past. The latter problem would require a theoretical model of the process. But, the proposed approach uses empirical models easily built from existing data and allows for the design of operating conditions for new products that lie close to past product grades. Literature Cited (1) Garcia-Munoz, S.; MacGregor, J. F.; Kourti, T.; Apruzzece, F.; and Champagne, M. Optimization of Batch Operating Policies. Part I. Handling Multiple Solutions. Ind. Eng. Chem. Res. 2006, 45, 7856–7866. (2) Garcia, S.; Kourti, T.; MacGregor, J. F. Model Predictive Monitoring for Batch Processes. Ind. Eng. Chem. Res. 2003, 43, 5929–5941. (3) Jaeckle, C. M.; MacGregor, J. F. Product Design Through Multivariate Statistical Analysis of Process Data. AIChE J. 1998, 44 (5), 1105– 1118. (4) Jaeckle, C. M.; MacGregor, J. F. Industrial applications of product design through the inversion of latent variable models. Chemom. Intell. Lab. Syst. 2000, 50, 199–210. (5) Garcia-Munoz, S.; Kourti, T.; MacGregor, J. F.; Matheos, A.; Murphy, G. Troubleshooting of an industrial batch process using multivariate methods. Ind. Eng. Chem. Res. 2003, 42 (15), 3592–3601. (6) Nomikos, P.; MacGregor, J. F. Multivariate SPC Charts for Monitoring Batch Processes. Technometrics 1995, 37 (1), 41–58. (7) Nomikos, P.; MacGregor, J. F. Monitoring Batch Processes using Multiway Principal Components. AIChE J. 1994, 40 (8), 1361–1375. (8) Eriksson, L.; Johansson, E., Kettaneh, N., Wold, S. Multi and MegaVariate Data Analysis: Principles and Applications; UMETRICS AB: Umea, 2001. (9) Kassidas, A.; MacGregor, J. F.; Taylor, P. Synchronization of batch trajectories using dynamic time warping. AIChE J. 1998, 44 (4), 864–875.
ReceiVed for reView October 23, 2007 ReVised manuscript receiVed February 26, 2008 Accepted February 27, 2008 IE071437J
