Ind. Eng. Chem. Res. 1992,31, 1665-1679
1665
PROCESS ENGINEERING AND DESIGN Robust Inferential Control for a Packed-Bed Reactor Hector M. Budman, Chris Webb, Tyler R. Holcomb, and Manfred Morari* Chemical Engineering 210-41, California Institute of Technology, Pasadena, California 91125 Two different inferential control schemes are applied to an experimental fixed bed methanation reactor. The first scheme, proposed initially by Brosilow, is designed based on Kalman filter estimation. The second less traditional design uses an estimator computed from the partial least squares regression method (PLS). The second approach was found to give superior performance when the nonlinear system under study is operated in a wide range of operating points. Due to the nonlinearity of the system it is essential to address the issue of robustness of the proposed schemes. This is formally done in this work using structured singular value theory. For the robustness analysis it is crucial to develop a realistic but not overly conservative uncertainty description. Since the PLS estimator uses a large number of measurements, a robust design based on the uncertainty associated with each one of the measurements would be very conservative. To overcome this problem, a lumped uncertainty description is proposed which is identified directly from experiments. 1. Introduction In many industrial processes the variables to be controlled cannot be measured fast enough to achieve desired control objectives or, in some cases, cannot be measured at all. However, readily available process measurements such as temperatures, pressures, and flow rates contain much process information. Exploiting relationships between these secondary variables and the variables to be controlled (primary variables) holds the potential to make more demanding control objectives accessible. For instance inferring the primary variables from faster secondary measurements will improve the performance of the closed loop system. This investigation focuses on the application of robust inferential control schemes to a packed bed methanation reactor. In addition to improved control, inferential control schemes often incorporate the benefit that secondary measurements can be leas expensive and more reliable than the primary measurements. Moreover, inferential control can add reliability to the closed loop control system. This is the case when the closed loop system is designed to use both fast and slow measurements in a cascade fashion, such as was initially proposed by Luyben (1973) and elaborated upon by Lee and Morari (1991). The cascade scheme will continue to function even in case of failure of either one of the measuring devices. Most of the published work on inferential control uses K h a n filter theory for the design of the inferential scheme (Brosilow and Joseph, 1978; Kumar and Seinfeld, 1978; Harris et al., 1980; Jorrgensen et al., 1984). These studies assume that the system is to be operated in a small neighborhood of one preselected operating point. Therefore, it is poasible to assume that in this neighborhood the process can be represented by a linear model. Using this assumption an optimal Kalman filter type of design can be implemented. Little work has been done on the application of inferential control to highly nonlinear systems
* Author to whom correspondence should be addressed.
operated in a wide window of operating conditions. In the present work we will evaluate different inferential control strategies when applied to a highly nonlinear system such as the chemical reactor for the case where the system is designed to operate in a wide range of operating points. The design of an inferential scheme consists of two parts: the design of the estimator and the selection of a controller. In the present work two different approaches to this design problem were investigated and compared. The first approach, which follows the Kalman filter type of design, assumes that the system can be modeled by a linear superposition of the effects of disturbances and manipulated variables on the primary and the secondary variables. Since this type of scheme was initially proposed by Brosilow (Brosilow and Joseph, 19781, it will be referred to in this paper as “Brosilow-type inferential scheme”. The second approach is based on the method of partial least squares regression, hence, will be described in this work as the “PLSscheme”. In the PLS scheme, in contrast with Brosilow’s scheme, no explicit superposition of the effects of disturbances and manipulated variables on the primary and secondary variables is assumed. Based on the results obtained using these two approaches, two main issues will be addressed 1. Accuracy of the inferential estimator. The goal of the inferential estimator is to accurately estimate the controlled variables. It will be shown that, due to the high nonlinearity of the process under study, the PLS technique results in significant performance improvement compared to the Brosilow approach. 2. Robustness of the inferential scheme. Due to the nonlinearity of the process and the wide range of operating conditions considered, studying the robustness of the proposed schemes, i.e., their operability in the presence of model uncertainty, is essential. In this work the robustness issues are formally investigated in the framework of structured singular value (SSV) theory (Doyle, 1982). The SSV framework assumes the uncertainty to be represented by norm-bounded perturbations to the transfer functions which model the process. Previously, work in
QSSS-5S85/92/2631-1665~Q3.QQ/Q 0 1992 American Chemical Society
1666 Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 1.1
1
0.5(, 0
,
,
,
, 20
,
,
,
,
,
,
,
40
,
,
, 60
,
,
,
,
I
80
Time (min)
Figure 1. Time response of exit concentration to step disturbance in wall temperature.
the K h a n filter framework used additive Gaussian noise in the system differential equations to describe the uncertainty, a less realistic description of model uncertainty. To our knowledge this work is the first experimental application of SSV theory to the design of a control system for a pilot plant chemical process, a fixed bed methanation reactor. Using this theory, it will be shown that a key issue for the analysis is the selection of an appropriate representation of the uncertainty in the system. By appropriate representation we mean realistic but not overly conservative. In the present application, a packed bed methanation reactor is operated under conditions of partial conversion of reactants. This experimental system is described in full detail by Webb (1990). The problem is to control the maximal bed temperature and the exit concentration of the reactants by manipulating the recycle flow rate and the heater power at the inlet to the reactor. The variables have to be controlled in the presence of changes in the reactor wall temperature. The control objectives are to obtain stability and good performance over a range of different operating conditions. If sufficiently many temperature measurements are available along the reactor, the maximal bed temperature can be obtained with sufficient accuracy from interpolation of some of these measurements. On the other hand, in industrial practice it can be very problematic to measure the reactant concentration on line. First, analyzers are relatively expensive to operate and require periodic calibration and maintenance. In addition, analyzers can have relatively slow sampling rates. The slow sampling rate results mainly from two factors: the time required to purge the line connecting the reactor exit to the analyzer and the time required for chromatographic analysis. In our system the purge time was approximately 6 min, and for our specific methanation reaction the time required to separate the gas mixture into its components is of the order of 10 min. The combined effect of these two delays results in a net delay of approximately 15 min. To get a feel for the effect of this measurement delay on the control performance of the closed loop system, we look at the time response of the reactant exit concentration to a step disturbance in the reactor wall temperature. This response is shown in Figure 1. The time constant for this response is of the order of the gas chromatographic measurement, Le., 15 min. Therefore, it is expected that the measurement delay will have significant impact on the performance of the closed loop if the gas chromatographic measurement is used for control. Although available at a slow sampling
rate, gas chromatographic measurements will not be used in this work for control. These measurements will be used solely for off-line comparison with the estimates obtained from the inferential scheme. Another approach for inferential control of a nonlinear process is to use an adaptive estimator (Gulandoust et al., 1988; Tham et al., 1991). However, there seems to be a certain reluctance in industry to employ adaptive schemes. This paper is organized as follows. In section 2 the inferential control problem will be presented. In this context the issue of measurement selection, studied by several researchers, will be briefly reviewed. In section 3 the experimental system will be described. In section 4 the design of the Brosilow-type scheme for the reactor is presented. In this section we will apply to the system a general measurement selection methodology, proposed by Lee and Morari (1992), which enables us to select the sensors for control and at the same time to address the robustness of the closed loop. The alternative design using PLS estimation and the robustness issues related to this design are presented in section 5. In section 6, the experimental results using each one of the two approaches are presented and compared. Final conclusions are presented in section 7. 2. Inferential Control and the Measurement Selection Problem Inferential control uses the values of secondary variables, temperature measurements in this study, and the manipulated variables to infer the effects of unmeasured disturbances on the controlled variables. The central problem in the design of an inferential scheme is the construction of an estimator which will produce accurate estimates of the variables to be controlled. Two main issues have to be addressed for the design of an inferential estimator: 1. Type of estimator. Inferential control studies can be divided according to the type of estimator used into two main categories: estimators designed using optimal Kalman filter techniques and estimators based on general regression techniques, e.g., the partial least squares technique. These two categories will be discussed separately in this paper. 2. Measurement selection problem. In practice, due to physical constraints and economic considerations,there is a limited number of sensors that can be introduced to the system. Therefore, the first problem that researchers confront during the design phase of an inferential estimator is that of measurement selection. The type, number, and optimal location of measurements for obtaining the best closed loop performance must be determined. This problem of measurement selection will be discussed for each one of the estimation techniques mentioned above. 2.1. Kalman Filter Estimation. Several investigations were conducted on the issue of inferential control using Kalman filter techniques (Brosilow and Joseph, 1978; Kumar and Seinfeld, 1978; Harris et al., 1980; Jargensen et al., 1984). Special emphasis was given in these works to the issue of measurement selection. These investigations can be grouped into two main categories according to how model uncertainty is accounted for in the analysis. In the first approach (Kumar and Seinfeld, 1978; Harris et al., 1980; Jargensen et al., 1984), mass and energy balances are expressed through partial differential equations. These equations are discretized and linearized around a nominal operating condition. The uncertainty in the system, due to nonlinear effects and uncertainty in the physical parameters, is modeled as a zero mean Gaussian white noise added to the linearized differential equations.
Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 1667 The main drawbacks of this approach are as follows: 1. The model uncertainty is lumped into a white noise Gaussian process noise term. This noise term is not a realistic representation of the uncertainty of the system due to process nonlinearities. Moreover, this type of description does not permit one to formally address the robust stability and robust performance of the closed loop system. 2. It assumes a very restrictive class of control structures. It is not clear how to extend the design technique to a different class of controllers, e.g., decentralized PID controllers which are widely used in industrial practice. 3. Substantial modeling effort and computational power are required to solve the resulting optimization problem. The second approach to the measurement selection problem describes model uncertainty as norm-bounded perturbations to the transfer functions describing the system (Brosilow and Joseph, 1978; Lee and Morari, 1991). This is in contrast to the previous approach in which uncertainty was represented by the process noise. This norm-bounded uncertainty is relatively easy to quantify from experimental data (Webb et al., 1989) and permits one to incorporate in a clear manner the issues of robustness (Morari and Zafiriou, 1989). This type of uncertainty representation was initially introduced, in the context of measurement selection, by Brosilow (1978). Brosilow’s selection technique holds only for steady state and for a specific structure of the model uncertainty. Lee and Morari (1992) extended the analysis of Brosilow by proposing a technique which enables one to conduct the measurement selection for general structures of model uncertainty and for the dynamic case. The application of this technique to the reactor will be presented in this paper. 2.2. Regression Techniques-Partial Least Squares Method. As mentioned in the Introduction, the objective is to design an inferential scheme which will result in good performance in a range of different operating conditions. The system under study, the fixed bed reactor, is known to exhibit highly nonlinear behavior in the selected window of operation. In all the aforementioned studies the estimator is designed on the basis of data obtained from a linearization of the system around a preselected nominal operating condition. Therefore, inaccurate estimation is expected when the system is operated ‘far” from this nominal condition. Partial least squares (PLS) uses a different approach to deal with the nonlinearity of the process. PLS was originally proposed by analytical chemists for “multivariable calibration problems” (Hoskuldsson, 1988; Wold et al., 1984) and applied lately by Mejdell and Skogestad (1989) to a distillation column for estimation of composition. Via the PLS approach a linear relation between the primary and both the secondary and the manipulated variables over a wide range of operating conditions is established. This has two direct benefits. First, this technique requires weaker assumptions than the Kalman filter technique regarding the linearity of the process. Use of the Kalman filter method assumes that the effects of disturbances and manipulated variables on both the primary and secondary variables obey superposition. In the PLS approach these assumptions are not required; the disturbances do not appear explicitly. Second, since the estimator is based on a range of operating conditions, changing from one operating condition to another does not require the development of a new estimator; from the outset, the need to function under different conditions is addressed. A main drawback of the
“en‘
H-0
Figure 2. Schematic of the experimental system.
PLS technique is that in its current form only static estimation is performed and it is not obvious how to extend the method to deal with dynamic information. In this work the dynamics of the estimator is corrected by augmenting the static estimator with a lead-lag compensator. Using the PLS technique, the accuracy of the estimator tends to increase as the number of measurements is increased. According to Brosilow increasing the number of sensors may result in increasing sensitivity to model error. However, it will be shown that this increased sensitivity is an artifact of the assumed uncertainty structure in the system. As was already mentioned, there is always a limitation on the maximum number of sensors which can be installed on a system. The PLS technique by itself does not solve the problem of sensor location. A simple algorithm for optimal sensor location to be used with the PLS technique will be discussed in section 5. In the present work we adopt the policy that a fixed number of measurements is available at fixed positions along the axis of the reactor and we are not free to either increase their number or change their positions. Therefore the sensors to be used for estimation with either one of the estimators discussed above have to be selected from this set of fixed sensors. 3. The Experimental System The reaction studied is the methanation of COz:
C02 + 4Hz
-
CHI + HzO
(1)
This reaction is exothermic and irreversible under normal operating conditions. It is catalyzed by a commercial nickel catalyst with which no significant side reaction has been observed. The reactor shown schematicallyin Figure 2 consists of a single 2.5 cm X 60 cm long stainless steel tube filled with a mixture of catalyst and inert alumina. Running down the center of the reactor is a 3-mm-wide stainless steel tube which contains 12 thermocouples at different axial positions, and surrounding the reactor tube is a bath of boiling Dowtherm. This fluid maintains a constant temperature at the wall of the reactor tube. This temperature is varied by manipulating the pressure of the bath. The flow rates of reacting gases, COzand H2,as well as the inert gas, N2,are individually controlled using mass flow controllers. The temperature of the resulting mixed
1668 Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 Table I. Oueratinn Conditions for Control Experiments amount units variable 0-1.5 5 % vol exit COz concn 268-290 “C exit temp 0-70 % of 1800-Wcapacity inlet heating power 12-20 slpm recycle flow rate 3.5 % vol inlet COPconcn 0.44 slpm COz feed flow rate 1.96 slpm Hzfeed flow rate 10 slpm N2feed flow rate 253-257 OC Dowtherm temp
stream is regulated using heating tape. The concentration of the product gas is determined with a gas chromatograph. The product stream is sampled every 3.5 min, and the results of the analysis are stored in a personal computer where they can be easily referenced off-line. A compressor connected in series to the reactor enables recycle of the products stream. The control objective is to regulate the exit concentration and the maximal bed temperature over a wide range of operating conditions. These operating conditions are listed in Table I and include both partial and complete conversion of the reactant COP. For the partial conversion conditions, the maximum bed temperature occurs a t the reactor exit whereas for full conversion a hot spot forms inside the reactor. The manipulated Variables are the recycle flow rate and the power supplied to the tape heaters at the inlet of the reactor. Increasing the recycle flow rate increases the space velocity in the reactor and makes the reactor behave (at least kinetically) more like a continuous stirred tank reactor (Hill, 1977). The methanation reaction varies from first order for low C02concentration to zero order for high concentration of this reactant (Van Herwijnen et al., 1972). On the basis of this fact it can be shown (Hill, 1977)that the steady-state C02exit concentration increases and the maximum bed temperature decream for increasing recycle flow rates. For the power supplied to the heaters, the situation is reversed. Increases in the power tend to increase the maximum temperature and decrease the exit reactant concentration. We considered disturbances in the reactor wall temperature only. Physically, this temperature is equal to the temperature of the boiling Dowtherm fluid which surrounds the reactor. There are two reasons for changes in the Dowtherm temperature: the pressure changes inside the Dowtherm container or the Dowtherm fluid does not boil uniformly. As a result, spatial temperature variations can occur due to a sudden increase of the flow rate through the reactor. It will be shown that the steady-state bed temperatures and the exit concentration change significantlywith small changes in the wall temperature. 4. Brosilow’s Inferential Control 4.1. Description of the Scheme. The main idea behind this scheme is to design an estimator of the disturbances and to combine it will a compensator which eliminates the effect of these disturbances on the primary variables. Let the system be described by the following dynamic model: s = G,dd + Gsmm c = Gcdd
+ Gcmm
(2) where s = vector of secondary measurements, c = vector of controlled variables, m = vector of manipulated variables, d = vector of unmeasured disturbances, and Gij =
___
I -
I
Figure 3. Schematic of Brosilow’s inferential structure.
transfer function relating output vector i to input vector
i.
A schematic of Brosilow’s structure is shown in Figure 3. For clarity of the analysis this inferential scheme can be put into an equivalent IMC (inter@ model control (Morari and Zafiiiou, 1989))structure. Q and E represent the IMC controller and the estimator, respectively. Ugng the IMC design procedure, is split into two parts, Q = Q,F, where F is a simple robustness filter and Q,is either and H2 or H, optimal controller designed for the nominal conditions. If G,, h e a stable right inverse, Q, = G,,-l is selected; otherwise, Q, is an approximate inverse of Gc,. Brosilow (Brosilow and Joseph, 1978)suggests for the case where dim s L dim d and for steady state the use of a least squares optimal estimator:
E = Gcd(0)GsdT(o)(Gsd(o) GSdT(o))-l
(3)
As dynamic estimator, Brosilow proposed either a Kalman filter or alternatively a simpler suboptimal estimator which consists of a static gain obtained from (3) augmented by a lead-lag compensator. The zeros and poles of this compensator have to be selected to minimize the integral square estimation error. For the special case where the number of independent measurements s is equal to the number of disturbances d , a dynamic estimator for c is obtained from (2):
E = G,.dG,d-l
(4)
when Gsd is assumed to be stably invertible. 4.2. Measurement Selection Procedure. Using the inferential control framework, Lee and Morari (1992) propose a new methodology for measurement selection and controller design based on stmctwred singular value theory. This analysis allows for structured norm-bounded uncertainty in each one of the system transfer matrices as well as frequency domain performance specifications. The results obtained using this method are much more general than Brosilow’s results, which hold under specific assumptions on model uncertainty and for the steady-state case only. In the present work, we will apply this measurement selection technique to the experimental reactor problem. Figure 4a shows the block diagram under consideration. In this figure
I
r
A u=
1 Al*
1
A,*
J
where Ai* are the normalized uncertainties associated with the transfer matrices which model the process. Lee assumed the uncertainty to be complex and dependent on frequency. This type of uncertainty gives rise to a family
Ind. Eng. Chem. Res., Vol. 31,No. 7,1992 1669
.
a
b
Figure 4. Block structure for steady-state performance analysis for Broeilow’s inferential scheme.
of linear time-invariant plants. The disturbances d and the controlled variables c are normalized by the frequency-dependent transfer matrices Wdand W,,respectively. The performance weight W, is selected to reflect the relative importance of the errors in the controlled variables. For simplicity, in the present work we will limit the analysis to steady-state performance only; the objective will be to minimize the steady-state offset in the controlled variables. The measurement selection is based on the minimizetion of the worst possible c l a d loop steady-state error over the considered uncertainty set A”. The block diagram shown in Figure 4a can be put into the M-A form shown in Figure 4b. The quantity cp*shown in this fis;ure is a scalar which multiplies the performance weight W,. Accordingly, a new output from the system is defined as c ” = c’c,*. Mathematically, cp* is selected as the largest scalar for which the robust performance criterion: max max
?
Ilc’llz < -
lld112 is satisfed. Then, for this configuration the worst possible closed loop offset over the considered uncertainty set A,, is d€D &€A.
(7)
Figure 5. Brosilow’s inferential control scheme applied to an uncertain system.
A,,*, in the transfer function Gsmand additive uncertainty, Asd*, in the transfer function Gsd. Because the uncertainty for the steady state and for the dynamic situations are estimated using different approaches, they will be discussed separately below. 4.3.1. Steady-State Uncertainty Description. The uncertainty bounds are obtained in the following manner: 1. The elements in the steady-state gain matrices, G,, and Gsd,are determined experimentally for a variety of different operating conditions. 2. For each element in these matrices, the static gain is bounded by the measured extremes. 3. The nominal model is computed by averaging these two extremes. 4. A bound on the uncertainty is obtained by computing the difference between the nominal gain and either measured extreme. This approach is relatively easy to implement experimentally. Uncertainty in G,,. The uncertainty of G,, is modeled as follows: where
for the estimator and the controller presented above. Considering the M-A structure in Figure 4b the value of cp*can be computed, for each measurement set, from the following equation using the structured singular value ( p ) [Morari and Zdiriou, 19891:
1
(12)
J.“
The computed bound on the worst possible offset, l/cp*, is used in this work for measurement selection. This bound is not tight since it is derived from a sufficient condition for robust performance. Therefore, the measurement selection criterion based on this bound may be conservative. The measurement set which maximizes c * has the smallest potential steady-state error and shoulcfbethe first set considered for control design. If cp*< 1,then the worst possible closed loop performance may exceed the performance criterion, as expressed by W,, and the measurement set may be eliminated from further consideration. On the other hand, if cp* > 1, robust performance at steady state is guaranteed. 4.3. Uncertainty Description. In order to account for the model mismatch, uncertainty in each one of the transfer matrices shown in Figure 3 should be considered. In the present work we consider, for simplicity, only two types of uncertainty (see Figure 5): additive uncertainty,
$A,) 5 1
~ ( A z 5) 1
(15)
where R = recycle flow rate, H = heating power, and T, = temperature measurement used for inference. The gains GT and G, are identified for different operating points wittin the selected window of operation (see Table 11). These gains are identified by exciting the system with a small step in either of the two manipulated variables while holding the other variable constant. Two sets of experiments are run. 1. The first set of experiments is changes in heating power holding flow constant-steps from 0 to 30% in the heating power for recycle flow rates of 10,12,15,and 17 standard liters per minute (slpm). 2. The second set of experiments is changes in flow rate holding power supplied to heaters constant-steps from 10 to 12,12to 15,and 15 to 17 slpm for 15% power. The
1670 Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 Table 11. Average Steady-State Gain with Corresponding Uncertainty for Transfer Functions GT.R, GTM,and GT.Tcorreln between
'4
?I b
6 s
2 3 4 5 6 7 8 9 10 11 12 13 14 15
-3.00 -3.08 -2.75 -2.58 -2.625 -2.75 -2.75 -2.75 -3.125 -3.25 -3.375 -3.125 -3.00 -2.75
0.5 0.08 0.25 0.08 0.375 0.25 0.75 0.75 0.875 0.75 1.125 0.875 0.5 0.25
1.1 0.96 0.84 0.7 0.66 0.64 0.6 0.56 0.53 0.5 0.44 0.36 0.21 0.13
0.1 0.08 0.08 0.06 0.06 0.04 0.04 0.04 0.05 0.02 0.04 0.04 0.13 0.17
2.16 2.33 2.33 2.5 2.58 2.91 2.74 2.74 3.16 3.16 3.33 3.5 3.5 3.33
0.5 0.33 0.33 0.195 0.25 0.08 0.08 0.08 0.16 0.16 0.33 0.16 0.16 0.33
0 0 0 0 0 0 0 0
+ +
0
-
-
nominal gains, 8 , and ~ 8, H , with their associated uncertainty bounds, WTz, and W T H , are computed as explained above and are presented in Table 11. For some choices of thermocouples A1 and A2 were found to be correlated. For example, it is observed that for thermocouples 13,14, and 15, for decreasing recycle flow rates, GTzR increases while GTzH decreases indicating a negative correlation. When this correlation is neglected, the resulting control design will be more conservative when compared with the design for which this correlation is considered. Uncertainty in GEd.Gadis a single-input single-output transfer function relating the reactor wall temperature to the selected measurement, T,. In this case the uncertainty was measured by imposing steps of -2 "C in Dowtherm temperature for two different operating points: 15% power, 10 slpm recycle, and 15% power, 17 slpm recycle. With these two experiments the system is excited at the extremes of the selected window of operating conditions. As is done for G,,, the steady-state gain and the uncertainty associated with that gain are computed from the results of these experiments. These results are also listed in Table 11. 4.3.2. Dynamic UncertaintyDescription. The robust stability of the closed loop is dependent on the uncertainty in the transfer matrix G,, but is not affected by the uncertainty in Gsd. Therefore, to test the robust stability of the system, we have to identify bounds on the dynamic uncertainty ASm.The elements of G,, are approximated by fmt-order transfer functions with dead time. The time constants and time delays are determined at different operating points from the same step response experiments conducted for steady-state identification. These values are presented in Table IV. Laughlin et al. (1986) have developed a technique for translating variations in transfer function parameters such as time constants and delays into an approximately equivalent norm-bounded uncertainty representation: P ( p ( I + WA); u(A) I 11 (16) where P is the "true" uncertain plant, p is the nominal plant, and W is the uncertainty weight. Using this method we compute frequency-dependent uncertainty weights w,d and WE,. 4.4. Steadystate Measurement Selection-Results. In principle, all possible combinations of temperature measurements should be considered when selecting the measurement set for inference. In the present work we limit the search to a single "best measurement". Beyond
I
I O
~ 0
,
,
,
, 5
,
,
,
,
~
,
,
,
,
,
,
,
,
15
10
Thermocouple number Figure 6. Performance bound cp*.
the obvious motivation to keep the estimator as simple as possible, additional arguments are as follows: 1. Changes in reador wall temperature are the only type of disturbance considered in the experiments. If the sensor noise is negligible, it is possible to obtain accurate estimates of a single disturbance using one single measurement with the estimator given by (4). Intuitively, one will expect that for a nonlinear system such as the one under study an increasing number of measurements may improve the estimation accuracy. In the next section it will be shown that this is not the case for a Brosilow-type estimator which uses a large number of sensors applied to a nonlinear process. 2. We observed that as more measurements are used for estimation it becomes more difficult to satisfy the robust stability criterion. To illustrate this point, we investigated the case where two secondary variables are considered for estimation. We examined 12 different pairings, each one including the maximal bed temperature and one of the 12 measurements along the reactor. A necessary condition for robust stability is that the following inequality is satisfied at steady state: (17) where MI, is defined by Figure 4a. None of the candidate pairs considered for estimation satisfied this condition. Following these considerations,one single sensor is used to infer the effect of the external disturbance on both controlled variables: maximal bed temperature and exit concentration. We can select the thermocouple to be used in the inferential control scheme based on the information in Table 11. The largest magnitude of the disturbances considered in the experiments is f 2 "C; therefore the disturbance weight w d is P(Ml1) I 1
Wd = 2 (18) For disturbances in the considered set the largest expected variations in the controlled variables are f0.5% in concentration and f15 "C in the maximal temperature. On the basis of these ranges of variation the weight W, is selected to penalize equally the expected largest changes in the two controlled variables as
w p = [0
OV15
]
(19)
The values of the worst steady-state offset, given by the
,
,
,
Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 1671 Table IV. Uncertainty Bounds for Time Constants and Delays time constant (8) delay (8) 661-701 200-600 164-264 0
Table 111. Transfer Function Models Based on Nominal Owrating Conditions R H Tw 0 . 0 0 3 8 ~ ~ -0.00034~-'~ -0.0095K2 C O 1 - 0.97452-' 1 - 0.97772-' 1 - 0.9578Z' 0.4386~-~ -0.444 0.0186~-'~ Tmcu 1 - 0.8431f' 1 - 0.9190z-' 1 - 0.86122-' 0.0195~-'~ 0.552r2 -0.547 T13 1 - 0.8301~-' 1 - 0.94302-' 1 - 0.85282-'
scalar l/cp* for each one of the measurements in the reactor, are computed and shown in Figure 6. The measurements in Figure 6 are numbered beginning from the reactor inlet. It is observed from Figure 6 that l/cp* decreases as thermocouple 13 is approached. After this thermocouple, a drastic increase occurs. This increase can be explained by the strong nonlinearity associated with the formation of the hot spot at the exit of the reactor. As pointed out in the previous section, negative correlation was found between the errors in ,G and,G for thermocouples 13, 14, and 15, respectively. For comparison cp* is computed for thermocouple 13 when the correlation is taken into account (cp* = 2.4) and for the case where the correlation is not considered (cp* = 1.7). This implies that by neglecting the correlation one overestimates the worst steady-state error by 30%. 4.5. Design of the Dynamic Compensator and Estimator. After the best single sensor for inference is selected, it is possible to proceed to the final step of the design, i.e., the design of the compensator and the dynamic inferential estimator. The transfer matrix of the process, G,,, is identified from step experiments. It is assumed that G,, can be modeled as a fmborder transfer function with time delay. The results of this identification are presented in Table I11 for a sampling time of 40 s. The compensator Q is designed using the IMC factorization presented by Holt and Morari (1985). This factorization procedure leads to a dynamically decoupled system with the minimal possible delays in the controlled variables. The matrix G,, is factored into an all-pass part and a minimum-phase part. The all-pass matrix, G ;,, is a diagonal matrix containing only dead-time elements. Using this factorization
r
-13
,
1 (21)
Table V. Operating Conditions for Steady-State Exwriments operating reijcle heat condition flow power wall temp mead concn no. (SlDrn) (%) ("C) (% vol) 1 17 5 257 0.74 17 50 257 0.24 2 50 254 18 0.69 3 0.17 19 65 255 4 13 0.18 45 256 5 15 0.26 45 256 6 18 0.99 7 20 255 19 1.41 10 254 8 1.3 10 19 9 255 15 257 19 1.06 10 1.18 18 11 10 256 0.81 19 35 12 256 0.93 20 35 13 255 0.53 16 35 14 255 0.38 17 30 257 15 0.31 15 20 256 16 15 0.52 256 5 17 14 0.65 254 5 18
Table 111. Using these models the estimator is computed from (4). For thermocouple 13 the estimator is given by r
1
0.8513
e, = (@,I-' D=
=
11
1 - 0.97772-'
0.444 1 - 0.843l.z-'
0.0038 1 - 0.97452-'
1 - 0.8528z-' 1 - 0.8612z-'
1
The robustness IMC filter has the form
F=
(1- e-TJX)z (2
- e-T./X)
I
By adjusting A, different controllers can be designed with different speeds of response. For the purpose of robust stability analysis, bounds for the time constants and delays of the elements of G,, are determined and presented in Table IV. These bounds can be translated into frequency-dependent uncertainty bounds using the already mentioned technique of Laughlin. The robust stability condition for the M structure shown in Figure 4b is given by ~ W 1 1 = )
1 - 0.9192-'
1
P(We,1QnJ'Ewe,J 5 1
V
w
(26)
From this inequality it is possible to compute the IMC filter time constant required for stability: ~ 2 8 3 s
I
0.0000706~~ (1 - 0.919z-'Xl - 0.9745z-') 0.00015 (1 - 0.843lz-'Xl- 0.97772-')
(23)
The zeros of D are inside the unit circle, and therefore, Qn is stable. The transfer functions G,, and Gsdare also modeled by first-order elements with time delays and are shown in
(27)
4.6. Discussion. Two main deficiencies were observed when using the previous scheme for the reador: inaccurate estimation and high sensitivity to model uncertainty. 4.6.1. Inaccurate Estimation. The objective of the estimator is to estimate the primary variables as accurately as possible, i.e., the exit reactant concentration and the maximal bed temperature. This objective has to be achieved for the nonlinear system when operated in a range of different operating conditions. To partially test if this objective is achieved when using the estimator proposed above, the reactor was operated in open loop around dif-
1672 Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 Table VI. Performance of Brosilow’s and PLS Estimators on Steadv-State Erueriments PLS Brosilow 1 12 1 12 m e a d operating sensor sensors sensor sensors concn condition no. (% vol) (% vol) ( % vol) (% vol) (% vol) 0.85 0.76 1 0.74 0.62 0.70 0.20 0.22 0.43 0.52 0.24 2 0.65 0.67 0.87 1.04 0.69 3 0.23 0.17 0.82 0.94 0.17 4 0.12 0.19 0.26 0.28 5 0.18 0.26 0.26 0.41 0.49 0.26 6 0.94 0.98 0.86 0.89 0.99 7 1.41 1.37 1.24 1.34 1.41 8 1.29 1.32 1.12 1.22 1.3 9 1.04 1.05 0.88 1.01 10 1.06 1.17 1.19 1.02 1.10 1.18 11 0.77 0.80 0.86 1.03 12 0.81 0.84 0.93 0.94 1.09 0.93 13 0.49 0.50 0.55 0.67 14 0.53 0.34 0.36 0.54 0.43 15 0.38 0.29 0.28 0.20 0.18 0.31 16 0.71 0.46 0.46 0.44 0.52 17 0.88 0.59 0.62 0.58 0.65 18 0.19 0.06 0.65 0.77 worst error 0.08 0.02 0.2 0.23 rms error
ferent operating points in the selected window of operation (seeTable V). At each one of the operating points defined in Table V, the temperature values obtained with the 12 thermocouples and the gas chromatographic readings of the concentration were recorded when the system reached steady state. We tested the accuracy of Brosilow’s estimator which uses sensor 13 for estimation. This measurement was found in the previous section to give the best steady-state performance. For simplicity we conducted the comparison for the exit reactant concentration only. Using the temperature values at thermocouple 13 and the correspondent values of the manipulated variables, the estimated concentration is computed from the linear equations of the system (2) from y = G,,6m + E(6T - G,,6m) + yn (28) where 6m = deviation in manipulated variables, 6T = temperature deviation at sensor 13, and yn = concentration at the nominal operating condition. Substituting (4) into (28) we find From the results shown in Table VI we observe that significant errors up to 0.65 90vol C02 occurred. We observe that the differences between the estimates and the actual values tend to increase as the system is operated “far” from the selected nominal condition which corresponds to 17 slpm recycle flow rate and 30% heating power. As mentioned in the Introduction, Brosilow concluded that for his choice of sensor noise and process noise, respectively, a larger number of measurements will result in more accurate estimation. It is interesting to test whether this conclusion still applies for the reactor operated at different operating points. For this purpose we compute a static estimator which uses all available 12 temperatures along the reactor according to (3). Table VI compares gas chromatographic readings for the different conditions with concentration estimates computed from (29). From the results we can conclude that overall the accuracy decreased, both the root-mean-square (rms)and worst errors increased as compared to the one measurement estimation case. However, certain im-
provements were observed for points “close”to the nominal operating point. Since Brosilow’s estimator was designed for linear systems with additive Gaussian process and sensor noise, the pattern of improvement “close” to one operating point coupled with performance deterioration “far away” should not be surprising. The effect of noise terms on the estimate can be expected to be filtered by this estimator as more measurements are used. However, as we demonstrated, with more measurements the estimator will not necessarily be more accurate when applied to a nonlinear system operated at different conditions. It has to be emphasized that (2) in effect assumes linear superposition of the effects of disturbances and manipulated variables on primary and secondary variables. This assumption can be arbitrarily inaccurate for nonlinear systems as demonstrated by the results in this section. This shortcoming is the motivation to look for an alternative estimator which will be more accurate for the nonlinear process studied in the present work. 4.6.2. Sensitivity to Uncertainty. In the Brosilow inferential scheme the controller is the inverse of the plant G,, at steady state. The condition number of G,, was found to be 812. This large condition number is due to the fact that the controlled variables are highly correlated at least in a part of the selected window of operation. This fact was found both experimentally in the present work and through numerical simulations in a previous work by Mandler (1987).In general this large condition number results in high sensitivityto model uncertainty (Morari and Zafhiou, 1989). This explains partially why the robust stability condition could not be satisfied even for a small number of measurements as was shown to subsection 4.4. An additional difficulty resulting from the high condition number is that the manipulated variables are easily driven to their saturation limits for relatively small external disturbances. For example, the heating power reaches its saturation limits for a change of the order of 3 deg in the wall temperature. Therefore an alternative dynamic compensator will be designed to avoid the aforementioned problems. This compensator will be used in conjunction with a partial least squares estimator. The design of the compensator and the PLS estimator are presented in the next section. 5. Partial Least Squares Estimation
In this section we will present the second scheme for inferential control for which the estimator is computed using the partial least squares regression technique. We will also address the issues of robustness related to this type of scheme. 5.1. PLS Estimator. The PLS design procedure seeks to minimize the squared error between measurements of primary variables and predictions of these primary variables made from secondary measurements for a set of different operating conditions. For this application, the primary variables were the outlet C02concentration and the maximal bed temperature. The secondary measurements were both manipulated variables and all 12 temperature measurements along the reactor. This estimator was based on the the first 14 operating points shown in Table VI; this data set is called the training set. One should note that only because steady-state data are considered, a static estimator is obtained. The remaining experiments in Table VI are used for testing the accuracy of the estimator. Knowledge of the process dynamics will be used later on to account for dynamic information. We found that a simple parabolic interpolation of the three highest contiguous temperature measurements along the reactor yielded good estimates of the maximal tem-
Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 1673 perature. Therefore, the PLS regression will be used solely to estimate the outlet concentration while the parabolic interpretation will be used to estimate the temperature. In general we expect that the accuracy of the regression will improve for an increasing number of measurements; hence, all 12 fixed thermocouples will be used for estimation. In the PLS technique the regression is performed directly from manipulated variables and temperature measurements to the concentration. This is perhaps the main reason that the PLS estimator is expected to produce more accurate estimates than Brosilow's estimator. For Brosilow's and K h a n filter schemes the transfer functions G,,, Gsd, Gcd,and G,, are separately identified and then the concentration is computed from (29) or a variant of (29)in the case of the Kalman filter. Thus,the coefficients for bm and 6T in the estimator are assembled from algebraic operations of the separately identified transfer functions. This is in clear contrast with the PLS estimator where these coefficients are computed directly from the regression procedure. The improvement in estimation accuracy resulting from this procedure will be clearly illustrated in this section. Before progressing, a review of the PLS regression algorithm is in order. The PLS algorithm was developed by Wold et al. (1984). A tutorial has been provided by Geladi and Kowaleki (1986), while Hijskuldsson (1988) has done a thorough job of analyzing the mathematical aspects of the algorithm. Partial least squares is very similar to principal components regression, which is in turn directly related to classical least squares linear regression. The regression technique is applied to the following matrix equation: y = Xr (30) In our case X is a matrix with the number of rows equal to the number of operating points considered for the regression. Each one of these row8 is composed of 12 temperature measurements and the value of the 2 manipulated variablea at each operating point as given in Table V. The vector y is composed of the measured concentration at each one of these points. The elements of the vector r are the weights which multiply each one of the temperatures and manipulated variables for estimating y. The classical least squares solution for r is r = (XTX)-lXTy (31) While this is the "optimal" solution, (31) often yields estimators that are very sensitive to noise if the condition number of the training set is high. A method to address this problem is principal components regression (PCR). PCR works by first using the singular value decomposition of X, X = UEV, to form a lower dimensional description of the input data. The columns of the square matrix V are eigenvectors of (Xrx). E is a matrix with the square root of the eigenvalues of (XTX) on the diagonal and zeros everywhere else. Moreover, E is organized such that each element uii L ujj for j > i. The elements of E are the square roots of sample variances of the data in the directions described by the columns of V. Directions with small eigenvalues have little variance and thus tend to have poor signal-to-noise ratios. Removing these directions is desirable because it tends to decrease the noise sensitivity of the predictor. To remove these directions, the rectangular matrix W is built from the first n columns of V,where the user must specify n. These columns correspond to the directions with the first n largest eigenvalues. A lower dimensional data set X,is obtained from the projection of the input data onto W, X,
= XW. Using this lower dimensional data set, the regression parameters are F = (WTXTXW)-'WTXTy
(32)
and the regression parameters for the full data are r = WF (33) The drawback of this approach is that the directions (W) are chosen on the basis of the correlation between the input variables (X), while the purpose of regression is to exploit correlations between the input variables and the output Cy). PLS addresses this directly by using directions of maximum correlation between X and y. Thus, PLS uses the same formula as above, but the elements of W are chosen differently than in PCR. Two weaknesses of both PLS and PCR are that one has to decide how many directions to use and that arbitrary scale transformations will alter the outcome of the regression. At present there is no theoretically general and accepted method for determining the number of directions to use. Most applications rely on enlightened trial and error approaches. Also, inputs with large variations will be favored over other inputs, regardless if this is reasonable. Thus, data must be properly scaled before applying these procedures. For the reactor data, 14 data points representing 14 different steady-state operating conditions were available. Each data point consists of 12 temperature measurements and 2 manipulated variables. Thus, the problem would appear to have an exact solution. However, the input data are highly correlated, making the problem very ill-conditioned. For instance, the matrix X in (30) has a condition number of more than 1OOO. To address this problem, the PLS approach was applied. First, the scaling issue needed to be addressed. A common suggestion in the literature is to scale all input variables to have unity variance. However, a significant source of error in the problem is noise in the thermocouples. Since the errors in each thermocouple are roughly the same, temperatures with the large variances should be favored. Thus, the two manipulated variables were scaled to have unity variance, and all temperatures had the same scaling factor applied; this single scaling factor was chosen so that the average variance of the scaled temperatures was unity. Next the issue of how many directions to use was addressed. When using PLS, this determination is typically done via cross-validation. In cross-validation, only a portion of the available data is used, building a trial predictor using a predetermined number of directions. The withheld data is then used to evaluate this trial predictor. As the number of directions is varied, the performance of the trial predictor on the withheld data (the testing set) will vary. Once the "best" number of directions has been determined via cross validation, one then uses all of the data to build the final predictor. This technique could not be applied directly for the reactor, since withholding any data would make the problem underdetermined. Thus, one could not determine if the predictor would improve if more directions were used. However, cross validation can still be used to determine a lower bound. If a certain number of directions are required to described well a portion of the data, then at least that many directions are needed to describe all of the data. The 14 data samples were divided into a training set of 7 samples and a testing set of 7 samples. The data were selected so that the training and testing sets had similar convex hulls; that is, the testing set involved little extrapolation from the training set. In underdetermined problems, PLS can still outperform classical methods. Using seven directions corresponds exactly to the mini-
1674 Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992
mum norm solution; however, six directions were found to give the best performance for the selected testing set. Having established that at least six directions were needed, an upper bound was sought. This was done by performing a principal component analysis of the data. Examining the eigenvalues of XTX revealed that 99.999% of the variation in the input data can be accounted for by using seven directions. Little benefit would appear to be gained from including more than seven directions. On the basis of this observation and the cross-validation results, the use of seven directions was deemed appropriate. With the number of directions determined, the PLS estimator was built using the the full (14 sample) data set. First, the seven-column W matrix was computed; next, W was used with (32) and (33) to generate a linear predictor. Combining this predictor with the scaling and mean centering discuwd above produced the desired predictor. Using this predictor, the estimated concentration is given by 9 = rxT + b (34) where r = (0.04673264 0.00530626
-0.05651377 -0.00775171 -0.05008534 -0.05812076 0.0322829 0.07138604 0.03232871 -0.02528921 -0.05728047 -0.01288784 0.02346117 4.005102123 (35)
b = 21.9773808 (36) The first 12 elements of x in (35) correspond to the 12 thermocouples, counted from the reactor inlet to the reactor exit, and the last 2 elements correspond to the recycle flow rate and the heating power, respectively. Since the above estimator involves the multiplication and addition of values of order 10 to produce estimates of order 0.01, eight significant digits are needed. This need for high-precision numbers in the estimator is not indicative of a "lack of robustness"; by scaling the inputs and outputs in the same manner as was done for building the predictor, low-precision values for the estimator parameters work just as well. We chose to use high-precision numbers in the predictor to avoid performing the centering and scaling of the input data for each new prediction. Whether one chooses to scale the input or to use highprecision numbers is an important implementation detail, but is not an issue of theoretical import. The accuracy of this PLS estimator is tested once again using the steady-state data presented in Table VI. The results of the estimation for the different operating points are shown in the same table. The worst error between the PLS estimation and the actual concentration is 0.06 9% vol in C02concentration. This is a significant improvement in the accuracy as compared to the Brosilow-type estimator. We also tested the conjecture that the accuracy of the PLS estimator may increase with an increase in the number of measurements. For comparison we compute a PLS estimator which uses one single temperature for estimation and compare it to the estimates obtained with the 12 sensor based estimator. Sensor 13 was selected since this is the sensor selected for Brosilow's estimator. The estimator is tested for all the operating points in Table VI. The worst error obtained is 0.19% COzconcentration, i.e., 3 times larger than the worst error obtained when all 12 temperature measurements were used for estimation. The main drawback of the PLS estimator is that it is based on steady-state information. Therefore, it is expected to produce good estimates at steady state but will be inaccurate for the dynamic situation. The formal ex-
1 l
0
'
"
'
2000
1
"
'
~
4000
I
~
~
6000
time (sec)
Figure 7. Time response of concentration to step in recycle flow rate for concentration: (-) measured with GC; (- - -) estimated with PLS estimator based on steady-state data, (- - -) estimated with PLS estimator augmented by lead-lag compensator.
tension of the PLS technique to deal with dynamic information is an open field of research. To test the dynamic accuracy of the estimator, the time response of the actual exit concentration to a step in recycle flow rate is measured (see Figure 7). The time constant for this response is of the order of 1200 s. If the concentration computed with the PLS static estimator is examined for this experiment, a time constant of approximately 300 s is observed (see Figure 7). This clearly indicates that the estimator dynamics obtained from the steady-state data is inaccurate. To correct this situation we propose to augment the PLS estimator with a simple lead-lag compensator: (37) The zero and pole of this compensator r A and TB can be obtained by fitting them to the actual data for one or more step responses of the system around different operating points as was done for the steady-state situation. For simplicity the fitting procedure was conducted for the step response to recycle flow rate shown in Figure 7 around the nominal operating condition (17 slpm recycle-30% heating power). From this fitting procedure we obtained rA= 300 s and TB = 1400 s. The time response for the concentration using the estimator given by (37) is shown in Figure 7. The estimator given by (37) has structure identical to that of the suboptimal estimator proposed by Brosilow as an alternative to the relatively complex optimal Kalman filter. Brosilow applied his inferential control to a distillation process. For this process he compared the closed loop performance using the optimal Kalman filter to the performance obtained using a suboptimal estimator, finding minimal differences between the two approaches. 5.2. Uncertainty Description. The PLS estimator proposed in this section uses a large number of measurements for estimation. This raises the question of an appropriate uncertainty description. Brosilow's analysis assumes that uncertainty is associated with each one of the measurements. Using this uncertainty description in the PLS scheme would be extremely conservative for the analysis, since the uncertainty regions for all the measurements would be added together. The worst possible combination of all these individual uncertainties would then define the uncertainty region for the predictor. Since the PLS scheme assumes that correlations between the secondary variables exist, such an uncertainty description would be inherently conservative and would become more
~
'
I
Ind. Eng. Chem. Res., Vol. 31, No. 7,1992 1675 the aforementioned technique of Laughlin to an uncertainty description suitable for robustness analysis. The uncertainty structure in Figure 8 is defined by ,. + Actm* Gci, = Gcfm (39) where / *
\
a
Figure 8. Block structure for steady-state performance analysis for the PLS inferential scheme. Table VII. Uncertainty Bounds for Gains, Time Constants, and Delays for the Transfer Functions Relating Manipulated Variables to the Estimated Concentration Using PLS and the Interpolated Maximal Bed Temperature gain time constant ( 8 ) delay (a) 1800-2200 0-160 G YH -0.007 to -0.015 G YR 0.06-0.16 1193-1507 0 500-1000 625-750 GT-H 0.12-0.3 GT-R -3.0 to -2.75 245-255 0
conservative as more measurements are considered. Therefore, an alternative uncertainty description is proposed which lumps together the individual uncertainty bounds into a more compact structure. The key insight leading to the lumped uncertainty description is to realize that the secondary measurements are only a means to an end. Uncertainties in these variables are only important inasmuch as they affect the overall accuracy of the estimator. Thus, only the uncertainty of the estimator is considered; the identification of the individual uncertainties is completely bypassed. The first step is to identify a transfer function which relates the manipulated variables directly to the estimated controlled variables. To identify this transfer function, the PLS estimator for concentration and the parabolic interpolation algorithm for the maximal temperature were implemented on the experimental system. Subsequently linear transfer functions relating the estimated controlled variables to the 2 manipulated variables and the disturbance are directly identified from step experiments. The resulting block diagram is shown in Figure 8a. The actual controlled variables are represented by c and the estimated variables by c'. To simplify the analysis we assume Gcm= G, (38) The procedure for identifying the uncertainty elements and the type of experiments conducted for this identification are identical to those discussed in section 4.3. As before, it is assumed that the system can be approximated by a first-order system with a delay. The nominal models for G&, and the variation in the time constants and delays along the window of operation are identified from step responses around different operating points (see Table VII). This parameter variation can be translated using
(43) Ac*m = diag (AI, A2, As, Ad) $Ai) 5
1
Z
= 1, 4
(44) (45)
where I",= is the maximal bed temperature and Y is the exit concentration. The actual value of the maximal temperature is identical to its estimate since these quantities are obtained using the same parabolic interpolation with three temperature measurements. Therefore c differs from c' only in the value of the concentration. The estimated value is obtained with the PLS regression while the real value is measured using the gas chromatograph. For the purpose of the analysis, an additive noise term n is added to the real concentration to account for the difference between the estimated and the real concentration due to regression errors. The magnitude of this noise is selected to be approximately equal to the worst error between estimatea and real values according to the results in Table VI, i.e., &0.06% C02concentration. 5.3. Design of the Dynamic Compensator. Using the inverse process model as controller was shown in section 4 to cause stability problems. These problems are due to the large condition number of the process transfer function. In order to correct this situation, a pragmatic decision is made regarding the performance requirements for the control scheme. Accurate control of the concentration requires integral action. However, if one only wishes the temperature to be maintajned at a reasonable level to avoid damage to the catalyst, then offset can be tolerated in this controlled variable. The concentration is paired to the recycle flow, and the maximal temperature is paired to the heating power. Since the concentration reacts much faster to recycle flow rate than to heating power, the above pairing is expected to give better performance for control of concentration than the reverse pairing. Explicit conditions were developed for the pairing problem for the case where integral action is required for all the loops (Morari and Zafiriou, 1989). Since integral action is not required in our case for the maximal temperature, these conditions do not apply to our problem. A more formal approach to the pairing selection is left for future investigations. A simple diagonal controller is
1676 Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992
proposed with proportional-integral control for the concentration and only proportional control for the maximal bed temperature:
to select the optimal location of the available measurements. The sensor location and the controller parameters which optimize steady-state performance may then be selected to solve the following optimization problem: min
min [(I/c,*)~+ n2]1/2 C
sensor location
\O
K2/
The expected offset between measured and predicted concentration, expressed by the noise n, is due to the regression error and is not affected by the selection of the controller parameters. On the other hand, the controller parameters can be selected to minimize the expected offset in the maximal temperature. The procedure used in section 4 for measurement selection is used in this case for selection of the controller parameters. As before, for the analysis the controlled variables c and the disturbance d are multiplied by a performance weight W, and a disturbance weight wd,respectively. The expected range of variation for the maximal temperature was shown in subsection 4 to be f15 "C. Then in order to scale the maximal temperature to 1 and to obtain integral action in concentration, the performance weight W, is selected as
-
w p = (0
)
O1/15
(47)
For disturbances of f 2 "C the disturbance weight is Wd = 2. The block diagram with the weights is shown in Figure 8a. This block diagram is transformed into the M-A structure shown in Figure 8b. Using this structure the controller parameters can be computed from a two-step iterative procedure: Step 1. Select parameters K1,r1,and K2which satisfy the robust stability criterion p(Mll) < 1. Step 2. Using the selected parameters from step 1, compute cp* from (8). The quantity cp* was shown in subsection 4.2 to be inversely proportional to the worst possible steady-state error. The controller parameters are selected to minimize cp* by iterating between step 1and step 2. Initial guesses for K1and T~ were derived from IMC tuning rules for single-input single-output control (Morari and Zdiriou, 1989) of the process represented by G m Since we found that the robust stability and the robust performance are mainly dependent on the value of the gain K2, the proposed two-step procedure is carried out by iterating on the value of this gain only. The controller parameters obtained from this procedure are K1 12 T I = 1350 K2 = 1 (48) For the above parameters cp* = 1.5. On the basis of the range of variation of the maximal temperature and for the selected performance weight the computed cp*is equivalent to an offset of f10 "C in the maximal bed temperature. In the present work we compute the PLS estimator using all the measurements available along the reactor. We assumed that their location cannot be changed and their number cannot be increased. For this situation we concluded on the basis of the steady-state results of Subsection 5.1 that the estimator which uses all available measurements is most accurate. If it was desirable to change the location of the sensors, the measurement selection technique proposed by Lee and presented in section 4 can be applied to the PLS scheme
(49)
where n2 is the variance of the noise representing the offset between the measured and predicted output. 6. Closed Loop Experiments
The main goal of the closed loop experiments is to test the stability and performance of the two inferential control schemes used in this work. By comparing the theoretical predictions of stability and performance to the experimental results, it will be possible to verify the identified uncertainty of the process. The stability and the performance results are presented separately in the sequel for each one of the schemes. 6.1. Experiments for Brosilow's Scheme. The experiments described below are conducted using thermocouple 13 for inference. Thus, both controlled variables, maximal temperature, and exit concentration are inferred from sensor 13. Stability Results. It is clear from Figure 3 that the stability of the closed loop depends only on the uncertainty in GSm.Therefore, it is possible to validate the uncertainty As,,, by comparing the experimental to the theoretical stability results. The IMC time constant required for robust stability as computed in section 4 is X 1 82 a. In order to check this result in the experimental system, we studied the disturbance rejection problem for different operating conditions using different filter time constants: 58,82, and 140 a. These experiments are conducted in the following manner: The pressure in the Dowtherm container is increased, resulting in a slow increase in wall temperature. The controller is activated at a predefined operating point. The pressure inside Dowtherm container is released, causing a step disturbance in wall temperature. The manipulated and controlled variables are monitored and the stability of the closed loop is determined. The operating conditions for which the experiments are performed and the stability results are presented in Table VIII. As is apparent from experiments 3, 4, 5, and 6 in Table VIII, it is more difficult to stabilize the system for lower recycle flow rates. That is, for small recycle flow rates the system appears to be most sensitive to model mismatch. By increasing the filter time constant, we are able to maintain stability at low flow rates, but a t the expense of a significantly reduced speed of response (experiment 1). According to Table I the system uncertainty was identified for flow rates between 12 and 20 slpm. However, using the filter time constant computed from the robust stability condition (82 a), instability was obtained for an operating condition corresponding to 13 slpm recycle flow rate and 15% heating power (experiment 2). This is not surprising considering that during the closed loop experiment the recycle flow rate becomes less than the minimal value of 12 slpm for which the procese was identified. This situation is clearly illustrated in Figure 9. Therefore we conclude that 88 long the manipulatad variables during the experiment remain inside the preselected window of operation, the stability condition obtained through p analysis is verified. Performance Results. The time history of the controlled variables associated with Experiment 5 are pres-
Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 1677 Table VIII. Closed LOOPInferential Control Exwriments robustmagniness filter tude of time disturb- stability init operating results const ance expt no. conditions stable 140s -1 "C 1 15% power; 13 slpm recycle; 256 "Cwall temp 2 15% power; 13 slpm 82 s -1 "C unstable recycle; 256 "C wall temp -1 "C unstable 3 15% power; 13 slpm 58 s recycle; 256 "C wall temp 4 15% power; 15 slpm 58 s -1 "C stable recycle; 256 OC wall temp -2 "C stable 5 5% power; 17 slpm 58 8 recycle; 257 "C wall temp -2 "C stable 6 15% power; 19 slpm 58 s recycle; 256 "C wall temp
b
2000
6000
4000 Time (5%)
8000
1.4,
;*80p Figure 10. Response of controlled variables to step decrease in wall temperature (experiment 5), using Brosilow's inferential scheme: (-) set point; (- - -) expected open loop deviation. 290
2 285
,
-
_________.___________._.___.._.__
$275
B 270
.c
0
6
1000
2000
Time (sec)
3000
4000
i;'++,
4000
2000
time
6000
8000
(sec)
='
...........................................
$12
,
,
I
I
,
,
,
,
,
,
,
,
,
,
Ei 0 6 e $04
0 - e o 5 0
2000 'Oo0
3000
4000
Time (sec)
Figure 9. Response of the manipulated variables to step decrease in wall temperature (experiment 2).
ented in Figure 10. Also shown in this figure are the set points and the expected deviations for open loop conditions. It is apparent from these figures that the inferential control scheme rejects the disturbance in wall temperature although it leaves significant steady-state offsets in both controlled variables. In the present work, for simplicity, only the steady-state performance criterion given by (7) is compared to the experimental results. More specifically, we check for all the experiments:
where l/cp* = 0.42 for thermocouple 13. Since the inequality was satisfied for all the performed experiments, we conclude that the estimated steady-state uncertainty is not invalidated. 6.2. Experiments Using the PLS Inferential Scheme. In this case the concentration is based on the PLS regression using 12 sensors and the 2 manipulated variable values while the maximal bed temperature is estimated using parabolic interpolation of 3 measurements. Stability Results. In this case the stability of the closed loop system can be used to validate the uncertainty Ad,,, (see Figure 8). The stability of the systems using the parameters given by (48) is tested for the conditions corresponding to experiments 4,5, and 6 in Table VIII. The system was found to be stable for operating conditions 5
2000
4000 time (sec)
6000
8000
Figure 11. Response of controlled variables to step decrease in wall temperature (experiment 5),using PLS inferential scheme: (-) set point; (- - -) expected open loop deviation.
and 6 but unstable for operating condition 4. This is due once again to the fact that for operating point 4 the recycle flow rate becomes less than the minimal recycle flow rate selected for identification of the system. Therefore, as concluded before, as long as the manipulated variables remain in the preselected window of operation for which the system was identified, the stability condition obtained with p analysis is verified and consequently the uncertainty ACrm is validated. Performance Results. The time evolution of the controlled variables for experiment 5 in Table VI11 is shown in Figure 11. A steady-state offset of 0.05% vol C 0 2 is obtained for the concentration while an offset of approximately 3 "C is obtained for the maximal temperature. The steady-state offset in concentration is reduced by a factor of 4 as compared to the result obtained with Brosilow's estimator. This is in agreement with the results presented in section 5 (Table IV),where the PLS estimator was shown to be overall more accurate at steady state along the considered window of operation. The estimated concentration using PLS regression was compared for experiments 5 and 6 to the actual concentration measured during these experiments. The results are shown in Figure 12. The maximal deviation obtained between the two variables is approximately 0.06% vol COP This is a very good accuracy considering that the accuracy of our gas chromatographic system is approximately *0.05% vol COP. These results partially validates the
1678 Ind. Eng. Chem. Res., Vol. 31, No. 7, 1992 z107
. i o 8 7 0
,
,
,
2000
I
4000 time (sec)
"
'
.
1
'
6000
'
"
~
8000
a
8 .= 1
O 6 1
051
,
,
,
,
2000
1
4000
~
'
"
I
'
6000
'
'
'
l
8000
time (sec)
b
Figure 12. Comparison between exit concentration estimated by PLS estimator (- - -) to concentration measured by GC (-) for (a) experiment 5 and (b) experiment 6.
dynamic correction of the PLS estimator using lead-lag compensation. 7. Conclusions Two inferential control schemes are applied to an experimental methanation reactor. The performance obtained with these schemes is compared and conclusions can be drawn from this comparison. The PLS estimator is shown to be significantlymore accurate for estimating the actual concentration in a wide range of operating conditions. A main disadvantage of the PLS regression technique is that it is based strictly on steady-state information. In this work it was shown that it is possible to correct the estimator to account for dynamic data using a simple lead-lag compensator. The zero and pole of this compensator are determined by fitting the estimate to experimental data. For this example, a PLS-based estimator was shown to perform better than another linear method, the Brosilow estimator. The reason for better performance of the PLS-estimator is that the construction of the Brosilow estimator is broken into two distinct steps while the PLS estimator is built to directly address the key performance issue, good prediction. In the first step of Brosilow estimator construction, the transfer functions G,,, Gad,G,,, and Gcdare separately identified. In the second step, these pieces are algebraically assembled to form the estimator. This second step relies heavily on the assumption that the plant is linear. The PLS-based estimator strives to provide the best possible (in the PLS sense) linear prediction without assuming that the underlying plant is linear. Thus, the PLS approach is preferred when one is attempting to build a linear estimator for nonlinear processes. An additional difficulty with Brosilow's scheme arises from selecting the controller as the inverse of the process model. Due to the high condition number of the process under consideration, this type of control will be very sensitive to model uncertainty. In order to overcome this problem, a decentralized controller is selected and its parameters are tuned for robustness using p criteria. The use of this controller in conjunction with the PLS estimator results in good closed loop performance. It was shown that for Brosilow's inferential scheme a larger number of measurements does not necessarily result in better estimation of the concentration for different
operating points. This is in contrast with the PLS regression scheme for which the accuracy was shown to improve as the number of sensors is increased. Developing a realistic but not overly conservative uncertainty description is crucial for analpis and tuning. In the PLS scheme, a lumped uncertainty between the manipulated variables and the estimated controlled variables is identified directly. A design using the uncertainty associated with each one of the measurements separately would be very conservative. Identifying correlations between the uncertainty elementa can reduce the conservativenessof the design. This type of correlationwas found in the context of the Brosilow scheme between the elements of As*. On the other hand, the PLS scheme apparently accounted for all the correlations between the parameters and no correlations were found between the elements of the uncertainty A,,,. A measurement selection technique developed by Lee and Morari in the context of Brosilow's inferential scheme was applied to the experimental system. The best measurement for estimation was shown to be the one closest to the hot spot but never crossed by it. The possible application of this measurement selection technique for the optimal location of the sensors to be used with the PLS estimator was explained. Future research encouraged by this work is a formal extension of the partial least squares technique to account for dynamic information and the application of nonlinear regression techniques for estimation.
Acknowledgment H.M.B. acknowledges the financial support of the Rothschild and Bantrell Foundations. T.R.H. is a recipient of a National Science Foundation Graduate Fellowship. Acknowledgment is made to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research. This research was supported in part by the Caltech Consortium in Chemistry and Chemical Engineering. Founding members of the Consortium are E. I. du Pont de Nemours and Company, Inc., Eastman Kodak Company, Minnesota Mining and Manufacturing Company, and Shell Oil Company Foundation. Nomenclature b = bias term of the PLS estimator in (36) c = primary measurement/controlled variable c' = controlled variable normalized by performance weight
WP
c" = product cc
*
d = unmeasurec! disturbances d ' = unmeasured disturbances normalized by performance
weight Wd E = inferential estimator F = IMC filter Gij = transfer function from variable j to variable i G :j = transfer function at normal operating conditions. K i= proportional control gain in eq (46); i = 1, 2 m = manipulated variable M = interconnection matrix for structured singular value (SSV) analysis defined in Figure 4b Mij = ijth element of M n = noise to represent the offset between measured and predicted concentration Q, = optimal controller designed for nominal operating conditions Q = product of Q, and F to form IMC controller r = regression vector i: = regression vector for the lower dimensional data set X, s = secondary variables
Ind. Eng. Chem. Res. 1992,31, 1679-1694
T,= sampling time U = matrix derived from the SVD of X V = matrix derived from the SVD of X W = projection matrix of X W,+ = weights used for SSV analysis described in Figure 4b; 1=1,2
Wd = disturbance weight W, = performance weight X = input data, each row corresponds to one input sample; each input sample consists of 12 temperatures and 2 manipulated variables X, = lower dimensional data set obtained from the projection ofXonW y = measured concentration yn = nominal concentration defined in (28) 6 j = deviation for the variable j A" = uncertainty matrix Ai* = individual uncertainty elements A = IMC filter constant p = structured singular value (SSV) Z = matrix of the singular values of X T~ = controller integral reset time in (46)
Literature Cited Brosilow, C.; Joseph, B. Inferential control of processes. AZChE J. 1978,24,485-509. Doyle, J. Analysis of feedback systems with structured uncertainties. ZEEE Roc. Part D 1982,129,242-250. Geladi, P.; Kowalski, B. Partial least squares regression: A tutorial. Anal. Chim. Acta 1986,185,1-17. Gulandoust, M. T.; Morris, A. J.; Tham, M. T. Adaptive estimation algorithm for inferential control. Znd. Eng. Chem. Res. 1988,27, 1658. Harris, T.; MacGregor, J.; Wright, J. Optimal sensor location with an application to a packed bed tubular reactor. AZChE J. 1980, 26. Hill, C. Chemical Engineering Kinetics and Reactor Design; Wiley: New York, 1977. Holt, B.; Morari, M. Design of Resilient processing planta-V The effect of deadtime on dynamic resilience. Chem. Eng. Sci. 1985, 40, 1229-1237.
1679
Hoskuldsson, A. Partial least squares regression methods. J. Chemom. 1988,2,211-228. Jorgensen, S.; Goldschmidt, L.; Clement, K. A sensor-location procedure for chemical processes. Comput. Chem. Eng. 1984,8, 195-204. Kumar, S.; Seinfeld, J. Optimal location of measurements in tubular reactors. Chem. Eng. Sci. 1978,33,1507-1516. Laughlin, D.; Jordan, K.; Morari, M. Internal model control and process uncertainty: Mapping uncertainty regions for SISO controller design. Znt. J. Control 1986,44, 1675-1698. Lee, J.; Morari, M. Robust control of nonminimum-phase systems through the use of secondary measurements: Inferential and inferential cascade control. Automatica 1992,submitted for publication. Lee, J.; Morari, M. Robust measurement selection. Automatica 1991,27(3),519-527. Luyben, W. L. Parallel cascade control. Znd. Eng. Chem. Fundam. 1973,12 (41,463-467. Mandler, J. A. Robust Control System Design for a Fixed-Bed Catalytic Reactor. Ph.D. Thesis, California Institute of Technology, 1987. Mejdell, T.; Skogestad, S. Estimate of process outputs from multiple secondary measurements. Proc. Am. Control Conf. 1989, 2112-2121. Morari, M.; Zdiriou, E. Robust Process Control; Prentice Hall: Englewood Cliffs, NJ, 1989. Tham, M. T.; Montague, G. A.; Morris, A. J.; Lant, P. A. Soft-sensors for process estimation and inferential control. J. Process Control. 1991,l (3). Van Herwijnen, T.; Van Doesburg, H.; De Jong, W. Kinetics of the methanation of carbon monoxide and carbon dioxide on a nickel catalyst. J. Catal. 1972,28,391-402. Webb, C. Robust Control Strategies for a Fixed Bed Chemical Reactor. Ph.D. Thesis, California Institute of Technology, 1990. Webb, C.; Budman, H.; Morari, M. Identifying frequency domain uncertainty bounds for robust controller design-theory with application to a fixed-bed reactor. Proc. Am. Control Conf. 1989, 1528-1533. Wold, S.; Ruhe, A.; Wold, H.; Dunn, W. The collinearity problem in linear regression: The partial least squares approach to generalized inverses. SZAM J. Sci. Stat. Comput. 1984,5(3), 753-743.
Received for review March 24, 1992 Accepted April 13,1992
Separation System Synthesis: A Knowledge-Based Approach. 2. Gas/Vapor Mixtures Scott
D.Barnicki and James R.Fair*
Separations Research Program, Department of Chemical Engineering, The University of Texas at Austin, Austin, Teras 78712-1062
A description is given for a prototype knowledgebased expert system, the separation synthesis advisor (SSAD),for synthesis of separation sequences for gas/vapor mixtures. The core of the SSAD is
the separation synthesis hierarchy (SSH),a highly structured, taak-oriented framework for representing separation knowledge. The hierarchy, based on interviews and information from the literature, emulates the approach that an expert process engineer follows. In ita current implementation, the SSH is limited to the preliminary sequencing of multicomponent gas/vapor mixtures using the following separation methods: (1)physical absorption; (2) chemical absorption; (3) cryogenic distillation; (4) membrane permeation; (5) molecular sieve adsorption; (6) equilibrium-limited absorption. Several examples of practical industrial separation problems are included. Introduction This paper is the second of a series on the development of a prototype expert system for the syntheais of separation sequences for fluid mixtures; the system is called the separation synthesis advisor (SSAD).Part 1 concentrates on separation system synthesis for liquid mixtures (Barnicki and Fair, 1990). Part 2 focuses on the parallel
problem for gaslvapor mixtures. The SSAD is a preliminary process design tool. Ita purpose is to formulate a limited number of feasible separation systems for a given multicomponent mixture. Final comparisons and optimization must be carried out with the aid of a process simulator, as the SSAD currently does not have the capability to perform a detailed economic analysis.
0888-5885/92/2631-1679$03.00/00 1992 American Chemical Society
