Ind. Eng. Chem. Res. 2010, 49, 8035–8047
8035
On the Calculation of Operability Sets of Nonlinear High-Dimensional Processes Christos Georgakis* and Lin Li Systems Research Institute and Department of Chemical and Biological Engineering, Tufts UniVersity, Medford, Massachusetts 02155
The calculation of the steady state operability characteristics of nonlinear high-dimensional processes is a difficult task computationally as it involves the mapping of high-dimensional, constrained input sets through the nonlinear process model to a high-dimensional output set. The reverse mapping, from outputs to inputs, also of interest, is equally challenging. In this paper we propose the selection of a finite number of points as an adequately accurate representation of the overall input-output mapping of the detailed model. This approach is motivated by the well established design of experiments (DoE) [Montgomery, D. C. Design and Analysis of Experiments; Wiley: New York, 2005. Box, G. E. P.; Draper, N. R. Response Surfaces, Mixtures, and Ridge Analysis; Wiley: Hoboken, NJ, 2007. Box, G. E. P.; Hunter, J. S.; Hunter, W. G. Statistics for Experimenters: Design, InnoVation and DiscoVery, 2nd ed.; John Wiley & Sons, Inc.: Hoboken, NJ, 2004.] methodology which has been successful in experimentally investigating similar relationships for processes that have not yet been modeled. For this reason, the proposed approach is called design of selective calculations (DoSC). The mapping of only selective points and the development of an interpolative response surface model for the output points enables the calculation of the desired operability sets with a much reduced number of calculations without any significant loss of accuracy. The applicability of the developed method is illustrated with two motivating examples and a plant-wide industrial process, the Tennessee Eastman challenge problem. Of particular interest is the demonstration that the response surface model with quadratic terms is rich enough to accurately describe even the complex phenomenon of input multiplicity. When input multiplicity exists, the image of the boundary of the input set does not fully describe the boundary of the image set. In such a case, mapping only the boundaries of the related sets is not sufficient. Besides its overall computational economy, the proposed method overcomes this important additional challenge. 1. Introduction 4,5
proposed a direct measure to Vinson and Georgakis quantify the operability of continuous processes. A process is called operable if the available set of inputs is capable of satisfying the desired steady-state and dynamic performance requirements, in the presence of the set of expected disturbances, without violating any constraints. The operability sets and the operability index (OI)4,5 are introduced as measures to quantify the process operability for both linear and nonlinear systems at the steady state. The related index can quantify the operability of existing plant designs by calculating the achievable output space [AOS ) M(AIS)]. Here the inputs take values in the available input set (AIS) and M represents the process model. This operability concept is also useful when designing a new plant as it can calculate the desired input set [DIS ) M-1(DOS)] that will enable the plant to achieve all the values in the desired output set (DOS). In the case that DIS is not a subset of AIS, this analysis can motivate changes in the initial process design to enlarge the AIS. Subramanian and Georgakis6,7 applied the operability framework to the analysis of reactors, units with strong nonlinear characteristics. They examined challenging lowdimensional nonlinear problems and they calculated the image sets AOS (or DIS) by tracing the boundary of the starting set AIS (or DOS). Parallel efforts in addressing the operability of continuous processes has also been undertaken by the research groups of Romagnoli8 and Lee.9,10 The primary objective of the present communication is to enable these operability calculations in high-dimensional nonlinear systems utilizing the steady state process model. For a * To whom correspondence should be addressed. E-mail:
[email protected].
case in which the input set is 2-D, one can quickly trace the 1-D boundary of the input set and calculate its image in the output set. In the case of a 3-D input set, the boundary is a 2-D surface and its image can only be calculated at a select set of points, say 11 on each edge, resulting in a calculation of 121 points in each of the 6 faces, for a total of 618 surface points. For higher dimensions, the number of such calculation points grows prohibitively large, requiring a new approach in defining the selective points at which to perform the calculation of the image points. Of interest here is to examine how many and which points will give us a reasonable approximation of the image set. The approach that we follow here aims at reducing the computational task. It is motivated by the design of experiments (DoE) methodology1-3 and uses the calculated response surface models to construct an interpolative but quite accurate model simplifying the exact process model. This makes the calculation of the operability sets for high-dimensional systems computationally feasible. For example, the calculation of the image of a three-dimensional cube can be obtained with just 33 ) 27 calculation points if we use the three-level design instead of the much larger number of points mentioned earlier. In the case we use the central composite design (CCD),1-3 the number of calculation points is reduced to 15. An optimal design11-13 can reduce the needed calculation even further. This approach has seen a significant number of applications in an engineering field other than chemical engineering where it is often referred by the term of “metamodeling”.14,15 The accuracy of the proposed method is demonstrated for several nonlinear processes, one of which exhibits input multiplicity, i.e., the map from inputs to outputs is not injective. We are interested in this characteristic because input multiplicity
10.1021/ie1009316 2010 American Chemical Society Published on Web 07/01/2010
8036
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
takes place in chemical processes and its occurrence presents a great challenge on the operability analysis. The challenge is caused by the fact that the boundaries of the available input set (AIS) do not fully map to the boundaries of the corresponding achievable output set (AOS). This makes the mapping of points only in the boundary of the starting set insufficient. Our approach resolves this challenge through the use of an estimated response surface model which approximates the exact model over the region where the calculations were made. The examples that we use to demonstrate the applicability of the proposed methodology are either square, with as many output as input variables, or with fewer outputs than inputs. However, there is no methodological limitation in addressing an additional important class of problems with a larger number of outputs than inputs. For such nonsquare systems, we previously introduced the concept of interval operability.16 In this case, one need only develop a few more approximate models, one for each output. This task does not represent any significant challenge beside the additional computational effort which is minimal. With such approximate models at hand, issues related to the process interval operability characteristics can be easily addressed. The structure of this paper is as follows. A summary of the prior work concerning process operability concepts is presented in Section 2. The proposed approach is detailed in Section 3. Two motivating examples are discussed in Section 4. In Section 5, a well-known plant-wide industrial problem with substantial complexity and high dimensionality is examined in order to test and verify the effectiveness of the proposed approach. This realistic industrial process exhibits input multiplicities which are accurately represented by the response surface model. Finally, we present our conclusions in Section 6.
input points is denoted as the desired input space (DIS), and we write DIS ) M-1(DOS). Finally, the expected disturbance set (EDS) represents the expected values of the process disturbances. Examples of these disturbances are process uncertainties such as operating parameters, uncontrolled compositions of feed streams, catalyst activity, or ambient conditions. With the operability sets introduced above, the servo operability index (OI) can be defined as follows:
2. Prior Work
In this section, we will introduce our approach of selective input-output calculations aiming to estimate an approximate response surface model. Detailed process models of complex processes often consist of heterogeneous model components including explicit equations, references to thermodynamic databases, as well as the use of legacy codes, among others. Such model character prevents one from performing necessary manipulations, such as the calculation of where the Jacobian is rank deficient. They also inhibit our ability to fully analyze the operability characteristics of the process. To overcome this difficulty, selective calculations will be designed and an approximate model of polynomial character will be estimated. In essence, we will be treating the detailed complex model in a parallel manner to the way we treat unmodeled processes where the method of design of experiments (DoE) is used to select at which conditions to perform the experiments and record the results. Here we also keep in mind that the cost of performing an additional calculation through the detailed model is much less than the cost of performing an experiment on a similar complex process. The design of experiments (DoE) methodology is a general strategy of planning, conducting, and analyzing experiments so as to obtain a good understanding of the response characteristic of the process for the minimum number of experimental runs. Such designs are based on sound statistical analysis, both for the design as well as the analysis of the results. An overview of the different designs will not be presented here as such information is widely available.1-3 Suffice it to say that, here, we will be initially interested in full factorial designs where each of the input variables (factors) will be assuming two (highlow) or, preferably, three (high-middle-low) values. The later type enables the estimation of models that include quadratic
In this section, a summary of the prior work on the concepts of steady-state operability is presented. As mentioned above, performing an operability analysis requires us to have a process model which relates the process inputs (manipulated and/or disturbances) to the process outputs. These process models can be described by the following generic form17 x˙ ) f(x, u, d) y ) g(x, u, d) h1(x˙, x, y, u˙, u, d) ) 0 h2(x˙, x, y, u˙, u, d) g 0
(1)
Here x ∈ Rnx are the state variables of the process, u ∈ Rnu are the control or manipulated input variables, d ∈ Rnd are the disturbance variables, and y ∈ Rny are the output variables. The model also includes equality and inequality constraints representing the process, product, and safety specifications, including the bounds on the magnitudes and the rate of change of the inputs. Here we will be concerned with a continuous process at steady state; therefore, the time derivatives in the process model will be set to zero. The operability sets defined by Vinson and Georgakis4,5 are summarized here. The manipulated inputs u are able to take values in the available input set (AIS). The corresponding set of values of the output variables y defines the achievable output set (AOS). We write that AOS ) M(AIS). Additionally, a desired output set (DOS) can be defined as the set of output points where we wish to be able to operate the process. The set of inputs required to reach the entire DOS can be calculated from the inverse of the process model. The collection of these
OI )
µ[AOS ∩ DOS] µ[DOS]
(2)
Here µ is a measure of the size of the corresponding set. In a two-dimensional space it represents the area and in three dimensions it represents the volume. The operability index indicates how much of the DOS is achievable with the available input set (AIS). The OI value is between 0 and 1. If the value of OI is equal to 1, the process is completely operable. If OI is less than 1, the initial expectations are greater than the process can achieve and some regions of the DOS cannot be realized. Besides the idealized reactor examples, Subramanian and Georgakis7 also examined the model of an industrial vinyl acetate reactor. This is a tubular reactor which exhibits input multiplicity. The occurrence of input multiplicity can be investigated by calculating the rank deficiency in the Jacobian matrix J, expressing the sensitivity of the output variables to the inputs.18 Besides the calculation of the images of the points in the AIS boundary, additional points in the interior of the AIS need to be imaged into the output set to fully calculate the boundary of the AOS. 3. Method of Selective Calculations
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
8037
terms. Additional designs of interest are the center composite designs in which one performs more selective calculations than in the two-level designs and much fewer than in the three-level designs, while retaining the advantage of being able to estimate quadratic terms in the response surface model. Depending on the design of selective calculations (DoSC) used, the number of calculations is fixed as well as the form of the response surface model and the maximum number of its parameters. In the case of two-level designs, the available data do not allow the estimation of quadratic terms and the response surface model (RSM) can only include linear and interaction terms (Note: For convenience, the symbol n here and in most of the following discussion is meant to imply nu.): n
y ) β0 +
∑
n
n
∑ ∑β uu
βiui +
(3)
ij i j
i)1
i)1 j>i
In the case of three-level designs of selective calculations, where high-middle-low values are given to each factor, the response surface model can include quadratic terms as well as third order interaction terms: n
y ) β0 +
∑βu
i i
n
+
i)1
n
n
∑ ∑β uu
ij i j
+
i)1 j>i
∑β u
2
ii i
ι)1 n
n
+ n
∑ ∑ ∑β
ijkuiujuk
i)1 jgi
(4)
k>j
Alternatively, one can use a second order RSM, which does not contain third order interaction terms and often provides a sufficient level of accuracy: n
y ) β0 +
∑βu
i i
i)1
n
+
n
∑ ∑β uu
ij i j
i)1 j>i
n
+
∑β u
2
ii i
(5)
ι)1
This is actually the most appropriate model for the center composite designs. In the response surface model in eq 4, one has 1 constant, n n linear terms (ui), ) n!/2!(n - 2)! interaction terms (uiuj), n 2 n quadratic terms (ui2) and 2 × third order term of the form 2 ui2uj. For n ) 3, 4, or 12, there are a total of 16, 27, or 157 terms, respectively, while for eq 5 one has 10, 15, or 91 terms, respectively. For the three-level full-factorial design, there are 33 ) 27, 34 ) 81, and 312 ) 531 441 selective calculations. One usually prefers to have a few more selective calculations of the detailed model than the numbers of parameters in the response surface model. Analysis of variance (ANOVA) calculations1-3 indicate which of the terms are insignificant and thus should be removed for the model. They can also offer a statistic about the fitness of the model to represent the data. To estimate the parameters of a quadratic response model represented by eqs 4 or 5, one can use the central composite design as it defines a smaller number of detailed calculation points. A 1/8 fractional center composite design involves 15, 24, and 537 selective calculations for 3, 4, and 12 factors, respectively. Of course, if we have 12 factors, the 537 selective calculations are too many for estimating the 91 terms of eq 5. In such a case, optimal designs are more appropriate (e.g., D-optimal).11,12,19
()
()
4. Motivating Examples In order to motivate the proposed approach to the highdimensional problem, we start with two examples of low
Figure 1. Comparison of AOS from the response surface model and the process model (nine points). The dashed line is the approximate AOS calculated from the RS model. The solid line is the exact AOS calculated from the process model.
dimensionality. The first is a single-CSTR process and the second is a system of two-CSTRs in series. In both cases, a first order reaction A f B is taking place. The two configurations have the same feed flow rates and final conversion specifications. We assume that both reactors are well mixed and operate open-loop without any other controllers besides the ones for level control. Both processes have been analyzed by Subramanian and Georgakis7,20 using the previous method of detailed mapping of the boundary points. The related equations for a single reactor are given in the Appendix. In order to make a valid comparison, the same operating conditions as well as the input and output variables are examined. For the operability analysis, the set-point of the volume hold-up of the reactor and the coolant flow rate are selected as input variables. The conversion and the reactor temperature are the output variables. The mathematical model is provided by Subramanian and Georgakis.7 Here we focus on a AIS region defined as the 50-100% range of the nominal volume hold-up and the 100-200% range of the nominal coolant flow rate. The full factorial design, with the two inputs as the factors and calculations at three levels, is selected to define the 9 () 32) points at which the input-output calculation will be made. The response surface model that is fitted to approximate the actual surface is given in eq 6. Here we select a quadratic model with the maximum possible number of quadratic terms. The model structure has 8 parameters to estimate for each of the two outputs from the available 9 calculation points. y1 ) R0 + R1u1 + R2u2 + R12u1u2 + R11u12 + R22u22 + R112u12u2 + R122u1u22 y2 ) β0 + β1u1 + β2u2 + β12u1u2 + β11u12 + β22u22 + β112u12u2 + β122u1u22
(6) The input variables ui in the above equation and in several subsequent equations represent the coded form of the input variables. The coded variables are a dimensionless transformation of the dimensional variable, input or output, so that the maximum and minimum values of the dimensional range correspond to the +1 and -1 values of the coded variables. From Figure 1, we see that the AOS calculated by the response surface model is very close to the one calculated by the point-by-point calculation through the detailed process
8038
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
Figure 2. Schematic of two CSTRs in series system with reaction A f B.7
model. This result demonstrates that the design of selective calculations holds the promise of yielding approximate models of noticeable accuracy in comparison to the detailed calculations, while offering a substantial ease in the analytic manipulation of the approximate model. This is an advantage but not a major one for the simple process examined in this example. However, the advantages of this methodology will be more substantial the more complex the process becomes and the more complicated the model is. The system of two CSTRs in series21 extends our approach to a slightly higher-dimensional problem. The schematic of this process is shown in Figure 2. In this case, we take the same ranges of the process variables as those of the single CSTR example. The difference is that we increase the number of inputs by introducing the coolant flow rate and volume of the second CSTR into the process. The inputs and output variables are defined as follows: (a) four (4) input variables, coolant flow rate, and volume of the first and the second CSTR; (b) two (2) output variables, conversion, and temperature in the second CSTR. As in the case of the single CSTR above, the range of the coolant flow rate is taken to be 100-200% of the nominal value, and the reacting volume in each of the two CSTRs ranges between 50% and 100% of the nominal value. For a 2-level 4-factor full factorial design, 24 ) 16 selective calculations are defined. The response surface model that can be estimated is given in the following equations. 4
y1 ) R0 +
∑Ru
i i
4
+
i)1 4
y2 ) β0 +
∑βu
i i
i)1
creation of a more accurate response surface model. The response surface model estimated is shown in eq 8 and the corresponding values of the parameters are given in Table 2. Here as well, the terms that are insignificant are shaded and their value is set to zero. Figure 4 depicts the approximate calculation of the AOS in comparison to the exact calculation. 4
y1 ) R 0 + 4
4
∑ ∑ R uu
ij i j
∑ ∑ β uu
(7)
ij i j
i)1 j)i+1
Here we need to estimate 11 coefficients for each equation from the available 16 data points, providing 5 () 16 - 11) degrees of freedom to assess the model adequacy. The comparison of the boundaries of the two AOSs which are calculated by the exact model and the approximate response surface model, eq 7, is depicted in Figure 3. The solid line denotes the AOS calculation from the detailed model and the dashed line the AOS calculation from the response surface model. The approximate model has only linear and interaction terms and is fitted at the vertices of the set. The estimated values of the parameters in eq 7 are given in Table 1. A p value of 0.05 is used as the margin of significance for a parameter. Those parameters that have a p value higher than the threshold are considered not significant and have been set to zero. Their entry in Table 1 is shaded. A more detailed design of selective calculations is considered, namely, the 3-level full factorial design in order to improve the accuracy of the AOS approximation. The 3-level 4-factor full factorial design has 34 ) 81 selective calculations of the detailed model. The large number of calculating points enables the
4
∑Ru
i i
4
+
i)1 4
∑ ∑ ∑R i>j
jgk k)1 4
y2 ) β0 +
ij i j
i>j j)1 4
4
4
∑ i)1 4
∑ ∑ ∑β
∑R u
2
ii i
i)1 4
∑∑
(8)
βijuiuj +
i>j j)1 4
ijkuiujuk
jgk k)1
+
4
βiui +
4
∑ ∑R uu +
ijkuiujuk
i>j
i)1 j)i+1 4 4
+
Figure 3. Comparison of AOS from the response surface model and the process model in 2-level 4-factor full factorial design. The dashed line is the approximate AOS calculated from the RS model. The solid line is the exact AOS calculated from the process model.
+
∑β u
2
ii i
i)1
In each of these two equations there are a maximum of 30 coefficients to estimate, initially assuming that all are significant, using the available 81 data points. Actually, one could easily substitute the above 3-level full factorial design with a fractional factorial design with 27 () 81/3) selective calculations. Alternatively, a central composite design defines 25 selective calculations. In both cases, one estimates a quadratic model with 21 terms depicted in eq 5. The values of the respective coefficients for the face-centered center composite designs are given in Table 3. We compare the relative precision of the two approximate models by their prediction of the temperature and conversion in the second reactor. The prediction point will be when all four of the coded variables are equal to 0.5. The full factorial design, based on 91 selective calculation points, predicts the temperature to be 73.986 ( 0.024 °C while the center composite design based on 25 selective calculations predicts 73.936 ( 0.14 °C. The corresponding predictions for the conversion in the second reactor are 98.505 ( 0.0078% and 98.496 ( 0.048%. The confidence intervals are based on a 95% confidence level. Both models are quite precise, with the center composite design giving a confidence interval that is 5 times larger. This is a very small price to pay for the benefit of having 56 () 81 25) fewer initial calculation points of the detailed model to be performed.
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
8039
Table 1. Values of the Significant Parameters in Equation 7 for the 2-Level Full Factorial Design
Table 2. Values of Significant Parameters in Equation 8 for the 3-Level Full Factorial Design index
0
1
2
3
4
12
13
14
23
24
34
Ri or Rij (°C) βi, βij, (%)
75.53 98.26
-0.10 -0.41
-1.39 0.51
-2.19 -0.06
-0.09 0.51
-0.30 0.10
0 -0.02
0.08 0.12
0.40 0.03
-0.02 -0.14
-0.09 0.01
index
11
22
33
44
123
124
134
234
112
113
114
Rii or Riij (°C) βij, βijk (%)
0.16 0.09
0.42 -0.14
0.73 0.02
0.01 -0.14
0.08 0.01
0 -0.03
-0.02 0
0.02 -0.01
0.07 -0.03
-0.04 0
-0.02 -0.03
index
122
133
144
223
224
233
244
334
344
Riij (°C) βijk (%)
0.08 -0.02
0 0.01
-0.02 -0.03
-0.12 -0.01
0 0.04
-0.13 -0.01
0 0.04
0.03 0
0.02 0
Table 3. Values of Significant Parameters in Equation 8 from the Center Composite Design index
0
1
2
3
4
12
13
14
23
24
Ri or Rij (°C) βi or βij (%)
75.53 98.26
0 -0.45
-1.44 0.51
-2.31 -0.07
-0.08 0.52
-0.31 0.10
0 0
0.08 0.13
0.39 0.03
0 -0.14
index
34
11
22
33
44
Rii or Riij (°C) βii or βiij (%)
-0.09 0
0 0.09
0.42 -0.14
0.70 0
0 -0.13
Another design of interest is the D-optimal design, an outgrowth of the work by Kiefer.11,22,23 This is a computergenerated design that selects the points where the experiments, and in our case the selective calculations, are performed so that a measure of the precision of the response surface model is maximized for a given number of experiments, calculations in our case. In the D-optimal design an algorithm selects the points at which the experiments (here the calculations) are performed so that the determinant of the information matrix is maximized. An additional reason that makes the D-optimal design, as well as other optimal designs, useful is the existence of constraints in the experimental design space or, in our case, the available
input set (AIS). A simple example of such a constraint on the input set could be that the sum of the cooling flow rates in the two reactors cannot exceed a maximum value. The number of experiments (or selective calculations) defined in each D-optimal design is determined by the response surface model that needs to be fitted. We have examined elsewhere24 two types of D-optimal design corresponding to two kinds of response surface models. The first model includes only the linear and the interaction terms between every pairing of the independent variables. The second model is a quadratic one. The number of regression coefficients to be estimated in these two response surface models are 11 and 15, respectively. Therefore, the number of selected calculations to be performed is at least 11 and 15, respectively. Having elucidated the proposed methodology with the above motivating examples, we now turn to addressing a more challenging industrial problem, the Tennessee Eastman process. 5. A Plant-Wide Industrial Problem
Figure 4. Comparison of AOS from the response surface model and the process model in 3-level 4-factor full factorial design. The dashed line is the approximate AOS calculated from the RS model. The solid line is the exact AOS calculated from the process model.
The Tennessee Eastman (TE) process is one of the most detailed dynamic simulations of a real process available in the open literature.25 The model of this process has been utilized as an example process in several publications in the area of process systems engineering.26-28 It describes the nonlinear transient material and energy balances in a five unit chemical process. There are 12 input variables and 42 output variables. Consequently, the available input set (AIS) is a subset in R12 and the achievable output set (AOS) is a subset of R42. Calculating the map of the AIS into the AOS is a substantial computational challenge. Assuming that the AIS is a 12-D orthogonal or oblique parallelepiped, the number of points that one would want to map without the benefit of the proposed design of selective calculations (DoSC) is 2 or 3 points on each
8040
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
of the 12 input dimensions, namely 212 () 4096) or 312 () 531 441) calculations, respectively. These are not difficult to perform, just time-consuming and, as we will show, unnecessary. What is much more difficult to do is to interpret the meaning of the 4096 or 531 441 42-D vectors that characterize the calculated points in the output space. Here are some issues that are not as easily resolved without the proposed DoSC methodology. (1) What is the minimum number and what is their preferred location of the input points to be mapped into the output space? Remark: The above-mentioned choice of two points in each dimension, 4096 total, is not appropriate in order to develop quadratic models. The choice of three points per dimension, resulting in 531 441 total points, is appropriate for quadratic model estimation. However, this many points are not necessary. (2) Are all the image points of the AIS boundary located in the boundary of the AOS? Alternatively, does the process have any input multiplicity characteristics? Remark: In such a case, segments of the boundary of AOS results from the mapping of points in the interior of AIS. These points are not easily identified without the proposed approach. Concerning the location of the calculation points in the AIS, let us, for a moment, assume that we will define them to be in the boundary of the 12-D parallelepiped. We can parametrize the boundary in a varying degree of detail. For this reason we first calculate that a 12-D parallelepiped has 212 () 4096) vertices, 12 × 211 () 24 576) 1-D edges, 66 × 210 () 67 584) 2-D faces, and so on until 12 × 2 () 24) 11-D boundary objects, serving a similar bounding role for the 12-D parallelepiped to that of the 6 faces for the 3-D cube (http://en.wikipedia.org/ wiki/Polyhedron#Vertex_figure). The minimum mapping task that we can do is for the 212 vertices. This however will not provide enough detail for understanding the AOS for processes that are nonlinear. For a more detailed input-output mapping calculation, we select p g 3 points along each edge. Then the number of calculations points is equal toN12,1,p ) 212 + 12 × 211(p - 2), where the first term represents the number of vertices and the second term represent the number points at the edges that are not vertices. If we parametrize in a similar manner, the 67 584 2-D faces and further on the 220 × 29 () 112 640) 3-D bounding cells (or cubes) we have the following number of input-output calculations, denoted, respectively, by N12,2,p and N12,3,p:
( ) ( )
12 12 2 + 12 12 12 N12,3,p ) 2 + 12 12 9 2 (p - 2)3 9
N12,2,p )
( )
( ) ( )
( ) ( )
12 11 12 10 2 (p - 2) + 2 (p - 2)2 11 10 12 11 12 10 2 (p - 2) + 2 (p - 2)2 + 11 10
(9)
Figure 5 depicts the dependence of the above two number of calculations, N12,2,p and N12,3,p, on the parameter p. In the same figure, we plot the value of p12, the number of the full factorial DoSC if p levels were to be considered for each one of the independent inputs. One can observe that the number of required calculations grows to be indeed very large, while the number of parameters of a given model remains constant. For a quadratic model it is 91 () 1 + 12 + 66 + 12). A systematic approach is needed in order to achieve the desired accuracy of the approximate model with a minimum number of calculations from the detailed model. In such a case the design of experiments (DoE) methodology can be very helpful in designing fractional designs of selected calculations (DoSC). Assuming
Figure 5. Dependence of the base 10 logarithm of the total number of calculation points on the number of points per edge.
that we wish to build a quadratic response surface model we need to have more calculations than the number of coefficients in the model. Such a model initially has one constant term, 12 linear terms, 66 interaction terms of the form xixj, and 12 quadratic terms of the form xj2, for a total number of 91 terms, as mentioned above. If we were to consider 12 independent inputs, the total number of calculations required for a 3-level full factorial design is 312 ) 531 441, a very tedious task to perform. Obviously 312 calculations are too many for the estimation of a set of 91 model parameters. If however, we perform a fractional design of selected calculations the task will be both more manageable to execute and the results much easier to interpret. Here we focus on the central composite design (CCD) that is always able to provide enough data for estimating a quadratic model. A 1/8 fractional CCD requires 536 calculations. This number drops to 281 in the 1/16 CCD case and to 105 in the case of a resolution V design1-3 with a minimum number of calculations or experiments. The point here is that an approximate quadratic model including a total of 91 terms, calculated above, can be estimated with as few as 105 in the above resolution V center composite design or 101 calculation points in a D-optimal design, instead of 4096 or 531 441. In our effort to develop approximate models, the DoE methodology greatly reduces the number of calculation points of the detailed model from which the approximate one will be developed. Not only does the DoE methodology drastically reduce the number of calculation points, it also determines in a systematic manner where these calculation points should be placed. Process Description and Prior Studies. The Tennessee Eastman process produces two products (G and H) from four reactants (A, C, D, and E), as well as an inert (B) and a byproduct (F) making the total components equal to eight. The reactions are as follows: A(g) + C(g) + D(g) f G(liq) G:product 1 A(g) + C(g) + E(g) f H(liq) H:product 2 A(g) + E(g) f F(liq) F:byproduct 3D(g) f 2F(liq)
(10)
All the reactions are irreversible and exothermic. The reaction rates are a function of temperature through an Arrhenius expression. The reaction to produce G has a higher activation energy than the reaction leading into product H. In addition, the reactions are approximately first-order with respect to each of the reactant concentrations. The process has five major
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
8041
Figure 6. Tennessee Eastman process with its control structure.33
processing units: reactor, condenser, vapor-liquid separator, recycle compressor, and stripper. Many control schemes have been applied to this problem.29-32 Lyman and Georgakis33 designed four closed-loop structures for the control consideration. Here we use structure no. 2 as our control structure because this structure is the preferred one of the four proposed structures and handles all the disturbances successfully. A schematic representation of the process and this control structure is depicted in Figure 6. Subramanian and Georgakis20 also performed an operability analysis on the TE process. They extended the operability framework to the analysis of plantwide systems. They defined the feasible operating region in the production-related variables as achievable production output space (APOS). The APOS of the process was calculated by solving an optimization problem which maximizes the production rate when the relative fraction of the two possible products is given. They also analyzed the effect of the key variables on the APOS limits. Their work mainly focused on observing the steady-state production capacity over a given range of interest in the product quality. Our work here concentrates on a more detailed operability calculation for the AOS. As we are interested in the steady-state operability, we set all the disturbances to zeros and run the simulation until the observed variables stop changing over time. Instead of the full problem of 12 inputs, only 4 are selected. They are the compositions of three major reactants (xA, xD, xE) in the feed of the reactor (stream 6) and the reactor temperature, T. This problem is not the largest size problem that can be handled.
Table 4. Base Values for the Input and Output Variable As Well As the Ranges of the Input Variable of the TE Process variables
base case value
range
units
A feed, xA D feed, xD E feed, xE reactor temperature, T
32.188 6.882 18.776 120.4
32-37 6.6-9.4 17-20 119-123
mol % mol % mol % °C
product G,xG product H, xH operating cost
53.724 43.828 170.6
mol % mol % $ h-1
However, it is complex enough to demonstrate the main point that we want to emphasize here. The product compositionsxG, xH in stream 11 and the total process operating cost are the considered outputs. Here the production rate is not changing, and we investigate the operability characteristics of the process and its control structure to produce different product compositions. Table 4 lists the base case values of the selected inputs and outputs. The reactor temperature and the reactant feed compositions xA, xD and xE have great impact on the product quality. The values of xG and xH are directly related to the product composition. We are also interested in the operating costs which are affected by the loss of raw materials in the purge stream. The TE process is constrained by the effectiveness of its controllers and the limits of the control variables. In defining the AIS, the following criteria need to be taken into consideration: (1) Within the low and high limits of the input variables, the TE simulation should operate without any shutdown. This indicates that no process violation occurs which can terminate
8042
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
the simulation. (2) The process should be able to reach a steady state after a proper period of time. This condition ensures the process to be stable in the AIS.
xG ) 60.24 + 7.76u2 - 0.73u3 + 0.74u4 + 0.089u1u3
All the calculations in this section correspond to the second control structure of Lyman and Georgakis.33 The range of the input variables, listed in Table 4, define the AIS selected, based on the above considerations. The feasible ranges of the input variables are relatively small due to the very sensitive nature of the process being open-loop unstable. We select the points at which we perform the calculation through the detailed model, following the procedure of design of experiments. We design a 3-level full factorial calculation for the 4 inputs, xA, xD, xE, and T. The design is defined through the values of the coded variables, where the 3-level full factorial design with 4-factors requires 34 ) 81 data points. The values of -1 and 1 in each calculation denote the lower and upper values of the corresponding input variable. The value of zero denotes the midpoint of the two bounds. For example, a calculation denoted by [-1, -1, -1, -1] corresponds to the following values of the four input variables [32, 6.6, 17, 119] and a calculation denoted by [1, 0, 0, -1] in the coded variables corresponds to the AIS dimensional point [37, 8.0, 18.5, 119].
-0.11u1u3 - 0.10u1u4 + 0.12u2u3 - 0.27u2u4 + 0.57u22 cost ) +176.70 + 10.54u2 + 7.75u3 + 2.42u4 + 0.85u1u2 +0.97u1u3 + 1.01u1u4 + 2.02u2u4 + 0.024u3u4 (13)
The above 81 simulations of the detailed model provide the corresponding values of xG, xH, and the operating cost. These data are used to calculate the three response surface models, one for each output, that aim to approximate the functional relationship between the input variables and the output variables of the detailed model. First a quadratic model is used as the response surface model. The mathematical form of the model is as follows: 4
y ) β0 +
∑ i)1
4
βiui +
∑
4
βiiui2 +
i)1
4
∑ ∑ β uu
ij i j
(11)
i)1 j)i+1
where the βs are the coefficients to be estimated and the us denote the input coded variables. The model is very precise resulting in an Radj2 value of 0.9998, 0.9998, and 0.9931 and in a PRESS value of 1.11, 0.84, and 91.38, respectively, for the RSM of the xG, xH, and the process operating cost. The corresponding quadratic models are given by eq 12, where all insignificant terms (p > 0.05) have been removed. xG ) 60.26 + 7.75u2 - 0.71u3 + 0.73u4 + 0.052u1u2 +0.075u1u3 + 0.049u1u4 - 0.16u2u3 + 0.27u2u4 - 0.62u22 +0.077u32 + 0.054u42 xH ) +37.37 - 0.038u1 - 7.69u2 + 0.63u3 - 0.67u4 -0.081u1u2 - 0.094u1u3 - 0.082u1u4 + 0.13u2u3 - 0.27u2u4 -0.036u3u4 - 0.045u12 + 0.60u22 - 0.079u32 - 0.052u42 cost ) +176.70 - 0.69u1 + 10.30u2 + 7.46u3 + 2.25u4 +0.59u1u2 + 0.69u1u3 + 0.78u1u4 - 0.47u2u3 + 1.74u2u4 -0.35u3u4 + 0.46u22 + 1.57u42
(12)
Now we will use a face-centered center composite design (CCD) that defines only 25 selective calculations from the detailed model. Almost one-quarter of the 81 calculations performed above. The corresponding RSMs are given in eq 13.
-0.13u2u3 + 0.26u2u4 - 0.58u22 xH ) +37.39 - 7.70u2 + 0.64u3 - 0.67u4 - 0.10u1u2
Here the model is almost identical to the previous one with the exception that some additional terms have been dropped as insignificant. However the dropped terms had a very small contribution in the previous model. This regression is characterized by the following values of Radj2 () 0.9996, 0.9995, and 0.9604) and of PRESS () 1.58, 1.81, and 132.27), respectively, for the RSM of the xG, xH, and the process operating cost. The Radj2 values are very good, though slightly smaller than before. Similarly the values of PRESS are very small, though slightly higher than before. The minor loss of accuracy in the three RSM models in this CCD design is amply compensated by the reduction in the needed detailed calculations from 81 to 25. We have timed the calculation of the detailed model at 16 points with the following coded coordinates ((0.5, (0.5, (0.5, (0.5) representative of the region we are concerned with. The execution time (All calculations with the detailed model and with the response surface model are performed using the GFORTRAN compiler of the GNU Fortran Project (http:// gcc.gnu.org/Fortran/) in a MacBook Pro (CPU, 2.4 GHz Intel Core 2 Duo, memory, 2 GB 667 MHz DDR2 SDRAM)) for all 16 calculations of the detailed model is 1.8582 min (111.490 s). Here each new steady state is calculated through the transient model of the TE process by integrating in time a set of ODEs, avoiding iterations of a steady state model. The execution time of the three response surface models in eq 13 is 0.98 × 10-6 s. (Because the execution time for evaluation of the three components of eq 13 is very small, 1 000 000 evaluations were performed to obtain a nonzero execution time.) This represent an execution time that is 1.14 × 108 times faster than that for the detailed model. AOS Calculation. In this section, we calculate the boundaries of the AOS corresponding to the boundaries of the AIS using the response surface model shown in eq 11. Note that we could as easily have used eqs 12 or 13 without a significant loss of accuracy. We also calculate the AOS via the detailed process model to obtain the exact AOS as a comparison. In order to be able to illustrate the AIS with three-dimensional plots, we fix one of the input variables at a specific value. The reactor temperature is chosen as the fixed input, at the values of 119 and 123 °C, which correspond to the -1 and +1 values of the coded valuable u4. The corresponding AOS plots are shown in Figure 7. The dashed lines represent the approximate AOS calculated via the response surface model, while the solid lines represent the exact AOS calculated by the detailed process model. Additional plots have been presented elsewhere.24 To understand in greater detail the characteristics of the AOS and to investigate the possibility of input multiplicities, we plot 2-D projections of the AOS by fixing two of the four inputs and letting the other two vary within their ranges. The outputs cost and xG are the selected two dimensions of the corresponding AOS, and the corresponding plots of the AOS are shown in Figures 8 and 9. The response surface model approximates the exact one quite well. In two additional figures,24 not presented here, the accuracy is similar. Most importantly, in Figure 9, and in one of the two not presented here, we also find that there exists intersection points in the image of the AIS boundary.
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
8043
case when there are three disturbances and the center composite design is used, there is a need for 15 selective calculations of the detailed simulation. In such calculations, the values of the control inputs uk are set at the center point values (uk ) 0). If we suppose that there are four controlled inputs, they require 25 points of detailed calculations for their RSM calculations. This brings the total number of selective calculations to 39 () 25 + 15 - 1); the center point is common in both RSM calculations. The overall model then has the form 4
y ) β0 +
4
∑βu
+
i i
i)1
4
4
∑ ∑β uu
ij i j
+
i)1 j>i
∑β u
3
+
2
iι i
3
∑∑
i i
+
i)1
ι)1
3
∑γd
3
∑γ d
γijdidj +
2
iι i
i)1 j>i
(14)
ι)1
This model has 24 () 1 + 4 + 6 + 4 + 3 + 3 + 3) parameters that can be easily estimated from the data of the 39 detailed calculations. In the second method, the disturbance variables are appended as additional variables to the vector of controlled inputs to define a longer input vector and the points of selective calculations is defined accordingly. For the example above with four controlled inputs and three disturbance inputs, we have a seven-dimensional appended input vector. If we denote with u5, u6, and u7 the three disturbance inputs, then the center composite design will lead into the following response surface model. 7
∑
y ) β0 +
7
βiui +
i)1
Figure 7. Comparison of exact AOS (solid line) vs approximate AOS (dashed line) at T ) 119 °C (a) and T ) 123 °C (b).
This implies that at least two different points in the boundary of the AIS lead to the same output point, the point of intersection, documenting the existence of input multiplicity. This has two important implications. (1) The simple quadratic response surface model proves itself able to represent such complex input multiplicity relationships of the detailed model. The occurrence of such model characteristics, as we will argue later, is difficult to discover in the detailed model. If the part of the detailed model consists of legacy code, then discovering the existence of input multiplicity might be close to impossible. (2) The calculated image of the AIS boundary is not sufficient to fully describe the boundary of the AOS, as it will be shown in the following. This has implications when we are interested to find out all the points of AOS, where we can operate the process. We have been mostly concerned here in building a response surface model between the controlled imputs and the process output. In the exact same way we can also build similar response surface models between the disturbance inputs and the process outputs. This can be done in two different ways. The first approach is to build a separate response surface model between the disturbances, di and the outputs yj. In the
7
∑∑
7
βijuiuj +
i)1 j>i
∑β u
2
iι i
(15)
ι)1
The advantage of the last model over the previous one is that interaction terms between the controlled inputs (u1-u4) and the disturbances (u5-u7) are estimated. A full size center composite design with 7 factors requires 143 selective calculations, a much larger number that the 39 for the first method. However the number of model parameters that need to be estimated are 36 () 1 + 7 + 21 + 7) so the 143 selective calculations are excessive for the fewer parameters to be estimated. This can be remedied in two different ways. The first one increases the complexity of the model by including third order interaction terms to the following form: 7
y ) β0 +
7
∑βu
i i
i)1
7
+
7
∑ ∑β uu
ij i j
j)i+1 i)1 7 7
∑ ∑∑
k)j+1 j)i+1 i)1
7
+
∑β u
2
7
7
ii i
+
i)1
βijuiujuk +
∑ ∑β u
2
ij i
uj (16)
j)i+1 i)1
This form has 99 constants to be estimated from the 143 points. Alternatively, instead of increasing the complexity of the model, fewer calculations can be defined by selecting a 1/2 fractional, a resolution V, or a small center composite design with 79, 45, or 37 selective calculations defined. Calculating Input Multiplicity through Response Surface Models. In the case of input multiplicity, the image of the AIS boundary is not sufficient to describe the boundary of the AOS. Then, one needs to calculate the images of a substantial number of interior points in the AIS, which will map into the boundary of the AOS. Not knowing a priori where the input multiplicity occurs, one needs to calculate the images of every point inside the AIS, a very tedious task in small dimensions and an almost impossible task in large dimensional problems. One needs a computationally effective way of accurately calculating all the AOS boundary points. We propose two methods, both using
8044
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
Figure 8. Exact AOS (solid line) and approximate AOS (dashed line) with fixed xE and T values and xA and xC varying in the intervals (32, 37) and (6.6, 9.4), respectively.
Figure 9. Exact AOS (solid line) and approximate AOS (dashed line) with fixed xA and xE values and xD and T varying in the intervals (6.6, 9.4) and (119, 123), respectively.
the estimated response surface model for an efficient calculation of the boundaries of the AOS. Method A. Having constructed the appropriate response surface models for the desired outputs y1,y2,. . .yny, as a function of the inputs u1,u2,. . .unu, we proceed as follows: (1) Delineate the image in the output space of the 1-D edge boundaries of the AIS using the response surface models. If additional details are needed, one can precede in calculation the image of the
2-D face boundaries of the AIS, and so on. (a) Remark: In the case that no input multiplicity is present, the above calculation is sufficient to define with sufficient accuracy the boundary of the AOS. (2) Fix p - 1 outputs except yq at values justified from the calculation in step 1. Then search for the maximum and the minimum values of yq using linear or nonlinear programming34-37 and note the corresponding input values of the uj’s. One of the optimizations might not be convex and
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
Figure 10. AOS calculation for TE problem with input multiplicity at xA ) 37, xE ) 20. The solid line represents the exact AOS. The dashed line represents the approximate AOS. The circles represent the maximal and minimal points.
global optimization methods might need to be employed.38 If the minimum and maximum values of yq are different from the ones calculated in step 2, the corresponding u point is an interior one in the AIS. (3) Repeat the calculation of step 2 for enough values of yq and rotate the selection of yq among the ny outputs. (4) Calculate the new boundary of the AOS by enlarging the calculations of step 1 using the additional points of steps 2 and 3 Here we focus on the first plot in Figure 9 in order to illustrate this approach. In this problem we have two inputs in the AIS, xD and T, which are varying in the intervals (6.6, 9.4) and (119, 123), respectively. We have fixed the other two input variables xA and xE as indicated at the top of each of the four subplots. In each subplot we plot the resultant values of the operating cost and xG. The exact calculation through the detailed model is given by the solid line and the dashed line is the result of the approximate response surface model calculated from a 3-level full factorial design of selected calculations. We observed that the shape of the AOS image in all four of the subplots of Figure 9 have a “crossover” point. This point is the image of two distinct points in the perimeter of the input set and implies the existence of an input multiplicity in the (xD, T) f (xG, operating cost) map. As a result there are interior points in the (xD, T) available input set that map into the boundary of the (xG, operating cost) image. We are applying method A above by fixing the value of xG and calculating the minimum and maximum possible values of the operating cost for all the available values of xD and T. In the framework of eq 12 or 13 we fix the value of xG and calculate the maximum and minimum values of the operating cost for all values of (u2, u4) in the (-1, 1) × (-1, 1) AIS, while (u1,u3) is equal to (-1, -1), (-1, 1), (1, -1), and (1, 1), corresponding to the four subplots of Figure 9. We are plotting one of these sets of calculations in Figure 10 for the case of (u1,u3) ) (-1, -1), where the circles denote the minimum and maximum values for 11 values of xG. These calculations achieve the full calculation of the AOS, in the case of input multiplicity using the response surface model. This calculation would have been much more difficult and likely impossible using the detailed model, depending on how convoluted the detailed model is. Although this method can calculate the true boundary of the AOS, the optimization might find a local optimum instead of a global one. Most importantly, this parametric calculation becomes a very difficult one as the
8045
Figure 11. Trajectory of det (J) ) 0 in the AIS at u1 ) u3 ) -1.
dimension of the problem increases. These limitations are addressed by method B. Method B. The input multiplicity occurs when the Jacobian matrix J, expressing the sensitivity of the output variables to the inputs, is no longer of full rank. This happens at the interior points of the AIS whose image becomes part of the boundary of the AOS.18 For the case of u1 ) u3 ) -1, the response surface model in eq 13 for xG and cost simplifies to xG ) 61.079 + 7.83u2 + 0.74u4 + 0.26u2u4 - 0.58u22 cost ) 169.22 + 9.65u2 + 1.386u4 + 2.02u2u4 (17) Then the Jacobian matrix is given by J(u2, u4) )
(
7.38 - 1.16u2 + 0.26u4 0.74 + 0.26u2 9.65 + 2.02u4 1.386 + 2.02u2
)
(18)
and det (J) ) 0 by det(J) ) 3.043 - 2.343u22 + 10.8u2 - 1.136u4 ) 0
(19) The trajectory in the (u2, u4) input space is shown as the line in Figure 11. Symbolic software can greatly facilitate the steps of method B, making it more attractive than method A. Method B is much more efficient than method A and holds great promise for high order systems. In such high dimensional problems, the critical question is how many dimensions are involved in the input multiplicity phenomenon. If the number of dimensions involved is small, method B can provide the answer. Otherwise an improvement might be needed. 6. Conclusions We have presented a new approach for the computation of the operability sets for nonlinear high-dimensional systems, denoted as design of selected calculations (DoSC). It is based on the methodology of the design of experiments and selects a relatively small number of points of calculations using the detailed model. An interpolative response surface model is then constructed, usually a quadratic one. This enables the efficient calculation of the necessary high-dimensional input-output maps.
8046
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
This approach falls under the concept of metamodeling used in other engineering contexts.14,15 The approach was demonstrated in two motivating examples consisting of a single CSTR as well as a pair of CSTRs in series. A more challenging case of a simulated industrial process was examined, the Tennessee Eastman challenge problem. The design of selected calculations used the full and fractional factorial designs, center composite designs, and optimal designs and in most cases a quadratic response surface model. In the processes examined, the quadratic response model was a sufficiently accurate approximation of the detailed model. Even in the case of the D-optimal designs where the number of detailed calculations is at a minimum, the accuracy and precision of the approximate model is sufficiently high. With such a response surface model in hand, the exact and approximate AOS were calculated. More interestingly, the response surface model approximated quite accurately the occurrence of input multiplicity. This enabled the accurate calculation of the response sets, including the calculation of the interior points in the input set that map onto the boundary of the output set. Two methods were introduced for this calculation. The second one, based in the calculation of the points at which the Jacobian matrix is not full rank, is the most efficient one. The examples that we use to demonstrate the applicability of the proposed methodology are either square, with as many outputs as input variables, or with fewer outputs than inputs. However, there is no methodological limitation in addressing an important class of problems with a larger number of outputs than inputs. For such nonsquare systems we have introduced the concept of interval operability.16 In this case one need only develop a few more approximate models, one for each output. This task does not represent any significant challenge beside the additional computational effort which is minimal. With such approximate models at hand, issues related to the process interval operability characteristics can be easily addressed. As the dimensionality of the problem increases, especially in the number of inputs, the proposed approach faces a strategic problem. This is related to the selection of the polynomial order of the response surface model but mostly on the calculation points needed from the detailed calculations, as well as their location. For a very large number of input variables, manipulated and disturbances, the center composite design is no longer an efficient design and one might have to rely much more on the Box-Behnken designs or one of the optimal designs, such as the D-optimal one used here. If, for example, we have 12 inputs, as the maximum number in the Tennessee Eastman simulation, the number of points the traditional CCD design needs is 4 289. In this case, the pure quadratic model has only 91 () 1 + 12 + 66 + 12) constants to estimate, which can be accurately estimated with fewer points. Even a 1/16 fractional center composite design defines 281 calculation points, which are 190 ( ) 281 - 91) more than the minimum needed. In such a case, one can select a minimum size resolution V CCD design which defines 105 calculation points, which are only 14 more than the number of coefficients to be estimated. In the case one opts for a more general model of cubic nature but without the pure cubic terms, the number of coefficients that need to be estimated are 377 () 1 + 12 + 66 + 12 + 220 + 66) and a 1/8 fraction CCD design is needed that defines 537 points for the location of the detailed calculations.
Table A-1. Parameters for the First-Order Reaction A f B parameters
values
parameters
values
β ξ γ φ δS δb
1.33 0.60 25.18 19.26 91.42 11.43
q0 x10 x20 x30 qc TR
1.00 1.00 -0.12 -0.12 6.56 599.67 °R
a liquid phase first-order, exothermic reaction of type A f B taking place in a jacketed CSTR. The mass and energy balances can be given in the dimensionless form, assuming constant physical properties and complete mixing as du2 ) q0 - q1 dτ d(u2x1) ) q0x10 - qx1 - φf(x2)u2x1 dτ d(u2x2) ) q0x10 - qx2 + βφf(x2)u2x1 - (δsu2 + δb)(x2 - x3) dτ dx3 ) qcµu1(x30 - x3) + ξµ(δsu2 + δb)(x2 - x3) dτ where x1 )
CA , CAR
T - TR , TR
x2 )
x3 )
TC - TR , TR
τ)
tQR , VR
A0e-γVR QR V (-∆H)C FCP R AR E γ) , µ) , β) , ξ) , RTR VC FCPTR FCCPC UAb UAS δb ) , δS ) FCPQR FCPQR γx2 f(x2) ) exp 1 + x2 q)
Q , QR
qC )
QCR , QR
φ)
( )
Here x1, x2, x3 are dimensionless variables corresponding to reactant A, reactor temperature, and jacket temperature, respectively. Since we are interested in the steady-state features, we set the time derivatives to zero. The resulting three algebraic equations can be combined to give a single equation representing the reactor energy balance with x2 the only unknown: g } q0(x20 - x2) +
βφf(x2)u2q0x10 q0 + φf(x2)u2 (δSu2 + δb)u1qc (x - x20) ) 0 qcu1 + ξ(δSu2 + δb) 2
Once all the relevant model parameters and feed conditions are specified, the above equations can be solved numerically for the steady-state temperature. Then from x2, it is straightforward to get x1 and x3 using the following relationships: x1 )
q0x10 q0 + φf(x2)u2
and x3 )
qcu1x30 + ξ(δSu2 + δb)x2 qcu1 + ξ(δSu2 + δb)
Table A-1 lists the parameters for the first-order reaction. Acknowledgment The authors gratefully acknowledge the financial support from PRF-ACS through Grant No. 45400-AC9.
Appendix: CSTR Model
Literature Cited
For the completeness of the present communication, we present here the equations that describe the CSTR model used. Consider
(1) Montgomery, D. C. Design and Analysis of Experiments; Wiley: New York, 2005.
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010 (2) Box, G. E. P.; Draper, N. R. Response Surfaces, Mixtures, and Ridge Analysis; Wiley: Hoboken, NJ, 2007. (3) Box, G. E. P.; Hunter, J. S.; Hunter, W. G. Statistics for Experimenters: Design, InnoVation and DiscoVery, 2nd ed.; John Wiley & Sons, Inc.: Hoboken, NJ, 2004. (4) Vinson, D. R.; Georgakis, C. A new measure of process output controllability. In Dynamics & Control of Process Systems 1998, Vol. 1 and 2; Georgakis, C., Ed.; Pergamon Press: Kidlington, Oxford, U.K., 1999; pp 663-672. (5) Vinson, D. R.; Georgakis, C. A new measure of process output controllability. J. Process Control 2000, 10 (2-3), 185–194. (6) Subramanian, S.; Georgakis, C. Steady-state operability characteristics of reactors. Comput. Chem. Eng. 2000, 24 (2-7), 1563–1568. (7) Subramanian, S.; Georgakis, C. Steady-state operability characteristics of idealized reactors. Chem. Eng. Sci. 2001, 56 (17), 5111–5130. (8) Bahri, P. A.; Bandoni, A.; Romagnoli, J. Operability assessment in chemical plants. Comput. Chem. Eng. 1996, 20, S787–S792. (9) Rojas, O. J.; Bao, J.; Lee, P. L. On dissipativity, passivity and dynamic operability of nonlinear processes. J. Process Control 2008, 18 (5), 515–526. (10) Rojas, O. J.; Bao, J.; Lee, P. L. A dynamic operability analysis approach for nonlinear processes. J. Process Control 2007, 17 (2), 157– 172. (11) Kiefer, J. Optimum Experimental Designs. J. R. Stat. Soc., Ser. B: Stat. Methodol. 1959, 21 (2), 272–319. (12) Atkinson, A.; Bogacka, B.; Zhigljavsky, A. Optimum Design 2000; Kluwer Academic Publishing: Dordrecht, The Netherlands, 2001. (13) Atkinson, A. C. Robust optimum designs for transformation of the responses in a multivariate chemical kinetic model. Technometrics 2005, 47 (4), 478–87. (14) Simpson, T. W.; Peplinski, J. D.; Koch, P. N.; Allen, J. K. Metamodels for computer-based engineering design: survey and recommendations. Eng. Comput. 2001, 17 (2), 129–150. (15) Sacks, J.; Schiller, S. B.; Welch, W. J. Designs for Computer Experiments. Technometrics 1989, 31 (1), 41–47. (16) Lima, F. V.; Georgakis, C. Design of output constraints for modelbased non-square controllers using interval operability. J. Process Control 2008, 18 (6), 610–620. (17) Georgakis, C.; Uzturk, D.; Subramanian, S.; Vinson, D. R. On the operability of continuous processes. Control Eng. Pract. 2003, 11 (8), 859– 869. (18) Koppel, L. B. Input Multiplicity in Non-linear Multivariable Control Systems. AIChE J. 1982, 28 (6), 935–945. (19) Atkinson, A. C. Examples of the use of an equivalence theorem in constructing optimum experimental designs for random-effects nonlinear regression models. J. Stat. Plann. Inference 2008, 138 (9), 2595–606. (20) Subramanian, S.; Georgakis, C. Methodology for the steady-state operability analysis of plantwide systems. Ind. Eng. Chem. Res. 2005, 44 (20), 7770–7786.
8047
(21) Subramanian, S.; Uzturk, D.; Georgakis, C. An optimization-based approach for the operability analysis of continuously stirred tank reactors. Ind. Eng. Chem. Res. 2001, 40 (20), 4238–4252. (22) Kiefer, J. Optimum Designs in Regression Problem - 2. Ann. Math. Stat. 1961, 32 (1), 298–325. (23) Kiefer, J.; Wolfowitz, J. Optimum Designs in Regression Problems. Ann. Math. Stat. 1959, 30 (2), 271–294. (24) Li, L. On the Calculation of Operability Sets of Nonlinear HighDimensional Processes. Ph.D. Thesis, Tufts University, Medford, MA, 2009. (25) Downs, J. J.; Vogel, E. F. A Plant-Wide Industrial Process Control Problem. Comput. Chem. Eng. 1993, 17 (3), 245–255. (26) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Integration of design and control for chemical processes: A review of the literature and some recent results. Annu. ReV. Control 2009, 33 (2), 158–171. (27) Bauer, M.; Craig, I. K. Economic assessment of advanced process control - A survey and framework. J. Process Control 2008, 18 (1), 2–18. (28) Maurya, M. R.; Rengaswamy, R.; Venkatasubramanian, V. Fault diagnosis using dynamic trend analysis: A review and recent developments. Eng. Appl. Artif. Intell. 2007, 20 (2), 133–146. (29) McAvoy, T. J.; Ye, N. Base Control for the Tennessee Eastman Problem. Comput. Chem. Eng. 1994, 18 (5), 383–413. (30) Ricker, N. L. Decentralized control of the Tennessee Eastman Challenge Process. J. Process Control 1996, 6 (4), 205–221. (31) Ricker, N. L.; Lee, J. H. Nonlinear Model-Predictive Control of the Tennessee Eastman Challenge Process. Comput. Chem. Eng. 1995, 19 (9), 961–981. (32) Sriniwas, G. R.; Arkun, Y. Control of the Tennessee Eastman process using input-output models. J. Process Control 1997, 7 (5), 387– 400. (33) Lyman, P. R.; Georgakis, C. Plant-Wide Control of the Tennessee Eastman Problem. Comput. Chem. Eng. 1995, 19 (3), 321–331. (34) Bazaraa, M. S.; Sherali, H. D.; Shetty, C. M. Nonlinear Programming: Theory and Algorithms, 3rd ed.; John Wiley & Sons: Hoboken, NJ, 2006. (35) Dantzig, B. Linear Programming; Springer-Verlag: New York, 1997. (36) Edgar, T. F.; Himmelblau, D. M.; Lasdon, L. S. Optimiztion of Chemical Processes, 2nd ed.; McGraw-Hill: New York, 2001. (37) Nash, S. G.; Sofer, A. Linear and Nonlinear Programming; McGraw-Hill Science, Engineering & Mathematics: New York, 1996. (38) Chachuat, B.; Singer, A. B.; Barton, P. I. Global methods for dynamic optimization and mixed-integer dynamic optimization. Ind. Eng. Chem. Res. 2006, 45 (25), 8373–8392.
ReceiVed for reView April 20, 2010 ReVised manuscript receiVed June 10, 2010 Accepted June 10, 2010 IE1009316
