Znd. Eng. Chem. Res. 1995,34, 1743-1754
1743
Identification of Uncertainty Bounds for Robust Control with Applications to a Fixed Bed Reactor Christopher J. Webb: Hector M. Budman,$and Manfred Morari* Chemical Engineering 210-41, California Institute of Technology, Pasadena, California 91125
A model-based robust controller is designed for a packed bed methanation reactor. To accomplish this objective, model uncertainty bounds are identified from experimental data. A physically motivated methodology of “regions mapping” was developed to compute the uncertainty bounds in the complex plane. This technique is compared to other existing nonparametric approaches for a simple nonlinear system and is shown to produce a more accurate description of the model uncertainty for the purpose of robust control design. This “regions-mapping” approach is then applied to a fmed bed reactor and uncertainty bounds are computed. A robust controller with a single adjustable parameter is designed for the reactor using internal model control (IMC) theory. The computed uncertainty bounds are experimentally validated using the IMC controller. 1. Introduction
The objective of this study was t o design a robust model-based linear controller for a packed bed reactor. The reactor is a highly nonlinear system due mainly to the nonlinear form of the reaction rate term. In general industrial practice, chemical reactors are operated and controlled over a variety of conditions. Therefore, if a linear model is selected to describe the reactor behavior, the parameters of this model will necessarily have to be varied to capture the nonlinearity and changes in operating conditions. A controller which is insensitive to this type of model uncertainty is defined as robust. Recent developments in control theory (Morari and Zafiriou, 1989) address the robustness problem in a formal mathematical framework. The uncertain process is modeled not by one single time-invariant model but by a family of linear time-invariant (LTI) systems. This is only an approximate representation of the uncertainty in view of the fact that one of the main sources of uncertainty in real processes are nonlinearities. On the other hand this approach allows the formulation of relatively simple but rigorous criteria to test the robustness of the control system. Mathematically, the family of LTI plants under consideration is parametrized by
where P ( w ) is the “true” plant, P(o)is a nominal LTI model, C is the multiplicative frequency dependent uncertainty bound, and A is any complex number with 1A1 < 1. While not strictly applicable to nonlinear or timevarying systems, linear robust control theory is often applied to systems with slowly varying or position dependent parameters. It is assumed that a controller which simultaneously stabilizes the plant for all possible parameter combinations will stabilize the true system. While it is sometimes possible to derive & for these systems from physical arguments, more often, the
* Author to whom correspondence should be addressed. Phone: (818)356-4186. Fax: 011-41-1-262-4362. E-mail:
[email protected]. t Current address: Icotron, Phoenix, AZ. t Current address: Chemical Engineering, University of Waterloo, ON, Canada. 0888-588519512634-1743$09.00/0
uncertainty bound is estimated empirically from observations of the system response. While much work has gone into developing robust control theory, relatively little is known about experimentally estimating uncertainty bounds for control purposes. This paper addresses this problem by proposing a technique for constructing the family of models .from input-output data and by showing the application of this technique to an experimental packed bed reactor. Currently, most uncertainty bounds are derived from bounds placed on the parameters of a parametric model (Postlethwaite et al., 1987; Palazoglu and Owens, 1987; Laughlin et al., 1986). Typically, a linear parametric model is proposed t o describe the nominal response of the system. The parameters used in the model are allowed to vary over specified ranges around their nominal values to account for deviations in the response of the true system from the nominal response. These parametric bounds can be determined either from experimentation or from physical arguments. When obtained experimentally, the system is generally excited at a number of different operating conditions and the dynamic response a t these operating conditions is compared to the nominal response. The parameters in the model are then allowed to vary in order to capture the different dynamics. While the parametric methods mentioned above are effective, the bounds obtained depend on the specific form of the system model. It is often desirable to have a black-box method which will directly identify a nonparametric dynamic model when given a set of inputoutput data and also give automatically a bound on the uncertainty associated with this model. Nonparametric models are of key importance in current chemical engineering practice due to the emergence of techniques such as DMC (dynamic matrix control) which uses this type of models (Cutler and Ramaker, 1980). This paper proposes a physically motivated methodology for identifying nonparametric models and their corresponding uncertainty bounds of single input-single output nonlinear processes such as the reactor under study. This method consists of two main steps: Step 1. By using spectral analysis, a series of nonparametric frequency domain models are identified from noisy input-output time records. These records are obtained from experiments conducted at different operating conditions. Step 2. By using the “regions-mapping” concept, a 0 1995 American Chemical Society
1744 Ind. Eng. Chem. Res., Vol. 34, No. 5, 1995
complex plane frequency dependent bound is found for all the individual frequency response models obtained in step 1. Spectral analysis techniques have been previously proposed for identifying uncertainty bounds (Loh et al., 1987; Rivera et al., 1987). Also the regions-mapping approach was used before for identifying uncertainty bounds for parametric models (Laughlin et al., 1986; Laughlin and Morari, 1989). The technique proposed in the present work is novel in the sense that it combines together concepts of spectral analysis and regions mapping to produce uncertainty bounds for nonparametric models. The bounds for nonparametric models are expected t o be less conservative than for their parametric counterparts. The paper is organized as follows. In section 2, spectral analysis theory for the identification of a nonparametric frequency domain model is reviewed. In section 2 our regions-mapping approach is presented. In the same section, two additional techniques based on spectral analysis and previously proposed in the literature for uncertainty estimation are also presented for subsequent comparison with our technique. In section 3, the estimates presented in the previous section are compared t o estimates based on Fourier analysis. It is shown that the uncertainty bound is actually composed of two parts: one part which measures the deviation of the model from the so-called “empiricaltransfer function estimate” (Ljung, 1987)and a second part which accounts for noise. In section 4, the nature of the input signal t o be used in the identification experiments is discussed. Two aspects of this signal are detailed: its frequency content and its shape. A simple technique is described for determining the frequency range over which the open-loop identification should be performed. This technique which uses a relay controller can, under certain circumstances, yield insight into the uncertainty bound. In section 5, the previously developed identification techniques are tested on a simple linear system in series with a nonlinear gain element. From this experiment, it is shown that the uncertainty bounds based on the residual spectrum can lead to misleading results when applied to a nonlinear system. The regions-mappingapproach overcomes some of the problems associated with the other methods and yields an uncertainty bound which is very close to the true bound. In the second half of the paper, the approaches presented in the first part are applied to the packed bed reactor process. This work is a first application of a robust controller, designed based on experimentally identified uncertainty bounds, to an experimental pilot plant reactor. In section 6, the control problem, the control of a hot spot in a fxed bed chemical reactor, is introduced. In this section, the flow rate through the reactor is chosen as the manipulated variable. In sections 7 and 8, the results of the identification experiments-closed-loop relay identification, section 7, and open-loop identification, section 8 are discussed. Finally, in section 9, a robust controller is designed and used to experimentally validate the accuracy of the uncertainty bound.
V
n
Figure 1. Block diagram of system with output noise.
bined with statistical theory to develop two different frequency dependent uncertainty bounds. These bounds differ only in their statistical assumptions. Then a new uncertainty bound is developed based on a regionsmapping approach. The dynamic experiments used to obtain this new bound are similar to the experiments mentioned in the Introduction for parametric modelling. Since spectral analysis theory is well established, the description which follows is only an overview of the theory as it applies to uncertainty estimation. The interested reader is referred to Ljung (1987), Brillinger (1975), Jenkins and Watts (19681, and Gardner (1988) for more details on the theory including the mathematical proofs. The system under consideration is shown in Figure 1 and is described by the convolution equation:
and its Fourier transform
where
This system is assumed to be stable, causal, linear, and time invariant. For simplicity, the signals entering the system are assumed to be bounded-magnitude boundedpower signals which may be either deterministic or stochastic. The manipulated variable u is a deterministic signal. The term u , which represents the effects of noise plus model uncertainty, contains both deterministic and stochastic components. In the analysis to follow, the expectation operator is implicitly applied to all calculations involving stochastic signals. If U,V, and Y are all known, then G can be deduced directly from eq 3. Unfortunately, V is usually not known and must be estimated from the identification procedure. Although Vis unknown, its effect on the outlet time series must be removed before the system can be identified. One way to remove the effects of V is to rely on correlation analysis. In correlation analysis, both sides of eq 2 are multiplied by u(t-2‘) and integrated over all time to yield
S_:(J:$(t-z)
U ( Z ) dr)u(t-T)
dt
+ S_:u(t)
u(t-T) dt
(5) 2. Spectral Analysis
In this section, the spectral analysis theory for the identification of a nonparametric frequency domain model is reviewed. It is shown how the residual spectrum obtained from spectral analysis can be com-
If u is uncorrelated with u , then and eq 5 reduces to
K,u(t) u(t-T)dt
A 0
Ind. Eng. Chem. Res., Vol. 34,No. 5, 1995 1745
e
problems arise. Given the series uN
= (u(Ts>,u(zTs>,...,u(NT,))
.YN
= (y(Ts)y(zT~),...y(NT,))
(where Ts is the sampling period) and their discrete then a simple estimate Fourier transforms UNand YN, of the transfer function at the points Wk = 2nk/NTs,k = 0 , (N/2)- 1is
(15)
Figure 2. Block diagram of system with output noise and uncertainty.
where (7)
The functions R,, and Ryuare called the autocovariance and cross covariance functions, respectively. [Strictly speaking, the covariance functions are only defined for stationary stochastic processes with mean value zero. The terminology used here is based on the terminology found in Ljung (Ljung, 19871.1 If the input is correlated with the disturbance signal, then computing the covariance functions effectively divides the disturbance signal into its correlated and uncorrelated parts. The correlated part then biases the estimate of g. For instance, if u ( t ) = J:,A(t - z) u(z) dz e(t) where e is a white noise signal (see Figure 2), then the estimated impulse response function is g A and not just g. The frequency domain equivalent to correlation analysis is spectral analysis. Taking the Fourier transforms of eqs 6-8 yields
+
+
Qyu(w)= G(w)Q u u ( w )
(9)
where
The overbar denotes the complex conjugate. Like the covariance functions, the functions Quu and Qyu are referred to as the autospectrum and cross spectrum, respectively. These equations play a central role in the identification of the plant, because
Using trigonometric interpolation for frequencies intermediate to the points W k , Ljung (Ljung, 1987) defines a continuous function, 6 ( w ) , which passes through the points 6 ( W k ) and calls it the empirical transfer function estimate, ETFE. For the given set of input-output data, 6 perfectly describes the data. By this it is meant that the output time series-can be reconstructed from the input time series and G without any error. However, due to its large sensitivity to noise, the ETFE can still be an inaccurate description of the process in the frequency domain; i.e., in the presence of noise the ETFE may still be a poor estimate of the true transfer function of the system. For most practical experiments the variance associated with this estimate is controlled by the signal to noise ratio and does not necessarily decrease as the length of the experiment increases (Ljung, 1985; 1987). To overcome these problems, 6 is smoothed. There are several ways in which this smoothing can be performed (see for example Kay and Marple, 1981) and Gardner (1988)). Two common techniques are Blackman-Tukey (“indirect”)smoothing (Daniell, 1946; Jenkins and Watts, 1968; Blackman and Tukey, 1959) and Cooley-Tukey (“direct”)smoothing (Bartlett, 1948; Koopmans, 1974; Hannan, 1970). In the Blackman-Tukey method, the spectrum (either cross or auto) at frequency wo is estimated using a weighted average of the unsmoothed spectrum around the frequency 00. In the Cooley-Tukey method, instead of averaging over frequency, smoothing is performed using a group averaging of several spectral estimates. The weighting function used for smoothing is termed spectral window in the frequency domain or lag window when the smoothing procedure is carried on in the time domain. For example, using the Blackman-Tukey procedure, the smoothed estimates of the transfer function and residual spectrum are
when QUu(w)f 0. The residual spectrum,-Qu,, is obtained by multiplying both sides of eq 3 by Y ( w ) Qu,(w>= @Jw)
- IG(w)12Q,,(w)
(13)
where 6$, and &E, are the smoothed estimates of the cross spectrum and autospectrum. These estimates are computed from 3
Unfortunately, when spectral analysis is applied to finite length discrete time series, certain covergence
N-1
1746 Ind. Eng. Chem. Res., Vol. 34,No. 5, 1995 1
N-1
reference books including Abramowitz and Stegun (1964). From these tables it is apparent that in the limit as a 0, f2,,,-2 00. Therefore, in order to be 100% confident that all uncertainty is bounded, the uncertain bound is infinitely large. In the present work we choose a = 0.05 which is equivalent to a 95% confidence that the uncertainty is being completely bounded. Clearly, for control purposes, it is desirable to describe the uncertainty with a 100% confidence. However, this cannot be allowed in a statistical sense since 100% confidence corresponds to an infinite uncertainty bound. Infinite bounds will clearly rule out the possibility of designing a realistic robust controller. A plausibility argument for using this bound in the design of a robust controller can be made by looking at the noise free linear system,
-
where w(Z) is the lag window and R(Z)is the discrete covariance estimate. There is always a trade-off between variance and bias: large spectral windows are needed for reduced variance while small windows are needed for reduced bias. In practice, a closing procedure is applied to the data whereby different size filters are used to smooth the data. The estimate which appears to be the best compromise between variance and bias is chosen as the smoothed estimate. Using the different smoothingtechniques, uncertainty bounds for the frequency response can be derived. Three different uncertainty bounds are presented. Two of them, previously reported in the literature, are strictly based on spectral analysis. The third one, developed in the present study, is based on a combination of spectral analysis and the regions-mapping concept. 2.1. Uncertainty Bound Based on Distribution Theory. Using the distribution theory associated with linear least squares estimations, Jenkins and Watts (Jenkins and Watts, 1968) convert the estimate of the residual spectrum into a bound on uncertainty. Specifically, these authors show that
1 - a (20) from which a confidence interval is derived
In these equations, Pr{x Iy> is the probability that x is less than or equal toy, f2,,,-2(1 - a)is the cutoff point associated with Fisher's F distribution of 2 and v - 2 degrees of freedom, and v is the number of degrees of freedom associated with the spectral window. Y is a function not only of the spectral window size but also of its shape. Thus the confidence intervals are a strong function of the amount of smoothing performed on the spectra. This is a key drawback for this computation since there is no straightforward method to decide on the right amount of smoothing. The empirical approach t o smoothing used in this work was window closing (Jenkins and Watts, 1968). This technique involves computing smoothed spectral estimates with a wide bandwidth window and then using progressively smaller bandwidth windows. This process is stopped when the bandwidth of the window is considered to be less than the smallest significant detail in the spectrum. Therefore, it is up to the designer and his physical knowledge of the system to decide which are the significant details of the spectrum. The F distribution is a sampling distribution which arises when computing the ratio of two independent x2 distributed random variables. Because it measures the ratio of two variance estimates, the F distribution is sometimes referred to as the variance-ratio distribution. Tables of the fvl,v2 are available in most mathematical
-
Y(u>= G(w)U ( w )
(22)
Let G+ be any transfer function estimate of G. Subtracting G+U from both sides of eq 22 yields
V(W)A Y(u>- G+(u)U ( U )= (G(u) - G+(u))U(U) (23) Multiplying by the complex conjugate
V(U)V ( U )= IG(u) - G f ( ~ ) l 2 U (U(U) u) (24) or
The similarity between eqs 21 and 25 becomes clear when [U(Y- 2)lf2,,,-2(1 - a)is set to 1. Using distribution theory, this factor predicts the fraction of the residual spectrum which is due t o model inaccuracy instead of noise. 2.2. Uncertainty Bound Based on Probability DistributionFunctions. A Werent uncertainty bound is obtained using a method presented by Loh, Corrga, and Postlethwaite (Loh et al., 1987). Using CooleyTukey smoothing, these authors obtain an asymptotically unbiased and consistent estimate of the spectrum. They then rely on the statistical theory of Wishart distributions (Goodman, 1963) to develop a set of probability distribution functions, PDFs, for the estimated squared coherency, pDyu12/@uu@yy;normalized transfer function gain; and phase. The coherency, I@ 12/@ @ can be viewed as the frequency depeniint corre ation coefficient between the input, u , and the output, y. These PDFs are implicit functions of the true squared coherence and phase. Once the PDFs are known, confidence intervals for both magnitude and phase can be derived in a straightforward manner. In this work the uncertainty bounds are approximated by circles a t each frequency in the Nyquist plane. In order to simplify the computations, the radii of the circles are selected equal to the computed magnitude bounds while the phase bounds are not explicitly used. The circle at each frequency obtained in this way is the one enclosed by the annular slice region defined by the magnitude and the phase confidence intervals at the specified frequency. Thus the uncertainty may be slightly underestimated since the corners of the annular slice are not included in the considered circular region. Because the developed probability distribution functions are rather complicated functions containing infi-
4 7 ,
Ind. Eng. Chem. Res., Vol. 34,No. 5, 1995 1747 nite sums of gamma functions, the cumulative probabilities are presented in graphical format. The confidence intervals are read directly from these graphs and depend on the amount of confidence one desires in the solution. Like the Jenkins and Watts bound, an uncertainty bound with probability of 1(100% certain) is infinitely large. As stated before, infinite bounds are of no use for robust control design. Hence, for this study, a 95% confidence level is assumed, and the uncertainty at this level fit to a power law relation for explicit evaluation. 2.3. Uncertainty Bound Based on Regions-Mapping Approach. An alternative technique for estimating the transfer function and its associated uncertainty bound can be derived from regions-mapping concepts. The first step of this technique consists of dividing the operating window, the admissible range over which y and u are allowed to vary, into a number of subsections. For each subsection, a local identification experiment is performed, and the resultant time series are analyzed using spectral analysis (eq 16). In this way, a set of local nonparametric frequency domain models is generated. At each frequency, these models represent a series of points on a Nyquist diagram. The second step of the technique is to find the smallest circle in the complex plane that bounds this set of points. In this way, a circular region is computed at each frequency. This computation step can be formulated as the min-max problem:
G= {local smoothed transfer function estimates}
G = {set of complex numbers) At each frequency, the transfer function estimate, em(wk), and the uncertainty bound, &(Ok)= max,+$ I@(wk) - & T w k ) l , are the center and the radius of the smallest circle which encompasses the set of plants &. An algorithm which solves eq 26 can be found in Laughlin and Morari (1989). In that work the minmax problem was solved to compute parametric uncertainty bounds, while in this work nonparametric models and their corresponding uncertainty bounds are identified. There are many similarities, and one important difference, between this technique and Cooley-Tukey spectral smoothing. In both cases, the input-output time series are divided into several smaller sections, a data tapering filter is used to smooth the data, and the transfer function for each section is determined. However, in the Cooley-Tukey method, the resultant transfer functions are averaged together to form a smoothed transfer function estimate, whereas with the regionsmapping approach, the local transfer functions are merely bounded and the transfer function estimate is taken as the center of the bound. While the regionsmapping uncertainty bound is not directly based on the residual spectrum, it is still a function of the amount of smoothing since the local models, @, are functions of the filter. If the uncertainty associated with the local models is relatively large, the uncertainty bounds for the individual estimates can be included in the global optimization (eq 26). As the number of local experiments increases, the likelihood of identifying the true bound, IG - 6-1,
increases. However, since the computed bound is a solution to a min-max problem, not all @ affect this bound. In practice, it is only necessary to identify the system at those points where the uncertainty bound will increase. While such a priori knowledge is not available in general, physical insight into the particular uncertainty present can be very useful for reducing the number of experiments. For instance, if the endpoints of the operating range are known to bound the frequency response, it is only necessary to identify the model at its endpoints. A big advantage of this scheme over the residualbased methods is the straightforward way in which it incorporates changes in the linear model (due to changes in the operating conditions)directly into the uncertainty description. These other methods must first relate changes in the model to an equivalent residual spectrum before these changes can be handled in the uncertainty description. As is demonstrated in section 5, these methods are not always effective at doing this. In addition, the residual lumps together both measurement noise and model uncertainty. While model uncertainty may cause a system to be unstable, noise does not. Therefore, a robust stability analysis based on the residual spectrum bounds may be inaccurate. With the regions-mapping approach, by splitting the range of operation into smaller ranges, the nonlinearity effects are approximately decoupled from the stochastic noise effects since in each one of the small ranges of operation the system is expected to behave in an approximately linear fashion. The traditional uncertainty estimates assume that the only model errors from an ideal linear plant are caused by stochastic noise. Therefore, these estimates will be more accurate when applied t o each one of the small ranges of operation rather than to the whole window of operation. An important question arising in the regions-mapping approach is related to the optimal number of subsections into which the admissible range of operation have to be divided in order t o obtain accurate uncertainty bounds. This topic is being left for a future investigation. This same philosophy of operating the system at several different operating points is also found in some of the parametric methods listed in the Introduction. However, unlike these other methods, the form for the parametric model which approximates the frequency response is not determined until after a bound on the uncertainty is found. This has the advantage of minimizing the size of & associated with the model. 3. Uncertainty Bounds Based on Fourier Analysis
In this section, the uncertainty bounds developed in the previous section are compared to bounds obtained from Fourier analysis. Consider again the system shown in Figure 2 for which the input-output time series, u,, and yn have been collected. Assume that 6, an estimate of G , is available and that the magnitudes of both E and A are bounded.
where E(wk) is the discrete Fourier transform of the noise term e and B E ( u k ) is a frequency dependent bound
1748 Ind. Eng. Chem. Res., Vol. 34,No. 5, 1995
on e. Then, one can ask the question, '"hat is the minimum bound on A necessary to describe the observed data?" This question can be mathematically posed as the minimization problem (29)
\O
otherwise
(33)
Using b as the uncertainty bound for the robust controller design, however, may prove to be too optimistic, for the noise signal is assumed t o perfectly line up with the residual and cancel its effects. On the other hand, an overly conservative bound on the uncertainty can be derived assuming that instead of cancelling the residual, the noise adds to it. In this case,
When using eq 33 (or 34),we are free to select the system model. In one extreme, G,can be taken t o be the ETFE in which case &(oh)A 0 independent of BE since by definition Y d w k ) - &wk) U d w k )= 0. However, the ETFE has already been shown to be a highly variable estimate of the true frequency response. In the other extreme, Ct may be a specific model. In this case, the model identification step has been completely avoided, and uncertainty may arise strictly due to a poor model choice. Somewhere between these two extremes are the models derived from spectral analysis. The estimated uncertainty bound is actually composed of two parts: one which measures the deviation of the modelAfromthe ETFE, IYdwk) - &ok) UN(wk)lI IUdwdI = IG - GI, and the other which measures the noise effects, B~(wk)/lUN(~k)l. When eq 33 is compared to eq 25, it is seen that the spectral techniques do not explicitly divide the residual between noise and uncertainty, but rather rely on statistics to implicitly make this factorization. This uncertainty analysis is actually a special case of a more general model validation analysis presented by Smith (Smith and Doyle, 1989). It has to be noticed that a main drawback for the computation of the bounds using Fourier analysis is that it is required to estimate the noise spectrum with good accuracy. For example, for the experimental system under study, the fxed bed reactor, it was not possible to estimate the noise accurately. 4. Input Design
This section addresses the design of the input signal used in the identification of an experimental system.
Two important characteristics of the input signal are discussed: its frequency content and its shape. 4.1. Frequency Content. The first step in any experimental identification is to determine the frequency range over which the plant is t o be identified. The input signal used to identify this system should concentrate its power in this frequency range. For control purposes, the frequencies of greatest importance are those around the crossover frequency of the system. While the crossover frequency can be determined from a number of different identification experiments, a closed-loop identification experiment using a relay controller predicts the crossover frequency directly and, potentially, yields insight into the magnitude of the uncertainty bound. A relay controller is a simple on-off controller for which the output of the controller is set to an upper limit, +d, if the input to the controller is greater than zero and to a lower limit, - d , if the input is less than zero. Using describing function methods (Vidyasagar, 1978;Astrom and Hagglund, 19841,it can be shown that the output of the plant in closed loop tends to oscillate at its crossover frequency. A simple variation on this controller is to precede the relay with a lead compensator (Smith, 1989). By adding phase lead to the plant, the crossover frequency of the plant with the lead compensator can be gradually increased. Eventually, a point is reached where too much lead has been added, and the output of the plant ceases to be periodic. There are two possible causes for this behavior. First, at higher frequencies the nonlinearities associated with the plant may dominate the response. For robust control design, this nonlinear response can be regarded as the output of a set of linear plants with complete phase uncertainty. Second, as the crossover frequency increases, the output of the plant becomes small. If the magnitude of the noise remains the same, eventually a point is reached where the relay is switching because of the noise and not the signal. We shall call the highest frequency for which the system under relay control is still periodic the breakdown frequency. As shown above, either noise or uncertainty is large at the breakdown frequency. From robust control theory it is known that because of stability and performance requirements, the closed-loop bandwidth is limited by the amount of uncertainty and noise in the system (Morari and Zafiriou, 1989). This limitation forces the closed-loopbandwidth t o be smaller than the breakdown frequency. For this reason, it is only necessary to identify the plant up to the breakdown frequency. 4.2. Shape. To reduce the amount of time needed in the identification experiment, it is desirable to have an input signal which simultaneously excites many different frequencies. One such signal is a pseudo random binary signal, PRBS (Godfrey, 1969). This signal is a two-level signal which uses a variable switching rate to excite multiple frequencies. When properly designed, the power spectral density of the PRBS is flat over the frequency range of interest and, therefore, approximates white noise in this range. However, since the amplitude of the signal is constant, the input can be carefully controlled to keep the system in a prespecified operating range. For the regions-mapping technique described above, the goal is to identify the system at a variety of different operating conditions including several different input
Ind. Eng. Chem. Res., Vol. 34, No. 5, 1995 1749 Table 1. Experiments Used To Identify Nonlinear Gain Model case
type of input
A B
PRBS PRBS PRBS PRBS at 3 different
C D
magnitude of input
method used to compute uncertainty distribution theory (section 2.1) probability distribution (section 2.2) regions mapping (section 2.3) regions mapping (section 2.3)
-3 to 3 -3 to 3 -3 to 3
PRBS -0.05 to 0.05superimposed
operating points
on input of fl,f0.333
levels. If a PRBS is used to identify the local dynamics for each operating condition, then the resultant input signal can be thought of as a PRBS superimposed on a series of different pulse functions where the size of the pulse determines the bias of the signal. Each pulse should be long enough to allow the PRBS to effectively identify the low frequency model response. The use of these two signals will be demonstrated in the next section where they will be applied t o the identification of a simple nonlinear system.
10
"True" Theoretical Upper Bound
5. Comparison of Uncertainty Bounds for a Simple Nonlinear System To test the techniques presented above, the uncertainty bounds of a simple nonlinear system were computed and compared to the actual bounds of the linearized model. This was done as a preliminary step before applying the same analysis to the reactor for which an exact nonlinear model is not available and then theoretical uncertainty bounds for comparison cannot be computed. It is shown that the bounds based on a single large PRBS inaccurately describe the true uncertainty and therefore may lead to misleading results. The system under consideration is a simple first order system in series with a nonlinear gain element. k ( t ) = -x(t)
+ u(t)
(35) (36)
For this system, the linearized transfer function between Y and U is given by
G(s) =
4
(2 - x0), s
1
+1
(37)
where the linearization is performed around x = X O . Confiningxa E [-1,11, the upper and lower bounds for the transfer function are given by
The gain of the identified nominal model is selected as the average of Glower and Gupper,
(39) with an associated uncertainty bound of
This system is simulated for 102.4 s a t a sampling rate of 10 Hz for a total of 1024 data points. Two different signals are used to identify this system. The first signal is a single PRBS with an amplitude of 3.
0.1
I
I
0.1
.1 Frequency (radlsec)
I 10
Figure 3. Comparison of upper bound on gain of transfer function for nonlinear system.
This amplitude is chosen so as to excite the state over the range x E [-1, 11. The second signal consists of a PRBS with amplitude 0.05 superimposed on a series of four equal length pulses with values of fl and 3~0.333. This system is identified using the various methods presented in Table 1. Figure 3 shows the predicted upper bound on the gain of the transfer function for all four test cases. In each case, the spectrum is filtered with a raised cosine lag window (Gardner, 1988). In experiments A-C, relatively poor estimates of the upper bound are obtained. All of these experiments use the single PRBS signal to identify the system, and therefore, it can be concluded that identification with this signal can lead to inaccurate results when used to determine the uncertainty bounds of nonlinear systems. Only experiment D, which uses the regions-mapping method, computes an uncertainty bound that is consistent with the true bound, Gupper. At low frequencies this bound deviates slightly from the true bound. This deviation is caused by a small amount of biasing that occurs in the local transfer function estimates. Figure 4 shows the calculated transfer function 8 , along with both uncertainty bounds. As expected, the lower bound fits the theoretical bound well, too. 6. Fixed Bed Methanation Reactor
The experimental model and uncertainty bounds identification of a packed bed reactor was conducted using the regions-mapping approach presented above. A schematic of this reactor is shown in Figure 5, and a complete description of the reactor can be found in Webb (1990). The reaction studied is the methanation of Con: CO,
+ 4H, - CH, + 2H,O
(41)
This reaction is exothermic and irreversible under
1750 Ind. Eng. Chem. Res., Vol. 34, No. 5 , 1995 10. 3
Frequency (radked)
Experimental Profiles
200
0 0.1
j.
10.
Frequency (radisec)
Figure 4. Frequency response of nonlinear model identified using regions-mapping technique. Superimposed on the frequency response are the theoretical and computed uncertainty bounds.
20
40 Distance from Reactor Inlet (cm)
Figure 6. Steady state axial temperature profiles for 3 slpm reactant. Table 2. Reactor Operating Conditions variable amount
coz Hz Nz inlet temperature Dowtherm temperature
Vent a
From
Heaters Dowtherm
Reservoir
Figure 5. Schematic of the experimental fixed bed methanation reactor.
normal operating conditions. It is catalyzed by a commercial nickel catalyst with which no significant side reaction has been observed. The reactor 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 used as an axial thermowell, 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 inside the bath. The flow rates of reacting gases, COz and Hz,as well as the inert gas, N2, are individually controlled using mass flow controllers. The temperature of the resulting mixed reactants is regulated using a heating tape. This temperature is held constant and roughly equal to the reactor wall temperature. The concentration of the
60
0.42 2.52 10-16 240 256
unitsa slpm slpm slpm "C "C
slpm, standard liters per minute.
product gas is determined using a gas chromatograph. The product stream is sampled every 20 min, and the results of the analysis are stored on a personal computer where they can be easily referenced offline. The operating conditions used in the reactor identification are presented in Table 2, and the steady state axial temperature profiles corresponding to these operating conditions are shown in Figure 6. The distinguishing feature in all of the profiles is the presence of the temperature maximum or hot spot. As the Nz flow rate increases, the height of the hot spot decreases, the width of the hot spot increases, and the position of the hot spot is shifted toward the exit of the reactor. All of these changes can be explained by the increased convection inside the reactor. This increase in gas convection tends to spread the generated heat over a wider area by transporting the energy away from the hot spot toward the exit of the reactor. For any given fixed operating conditions, the location and height of the hot spot are unique functions of the N2 flow rate. It is desirable to maintain the hot spot near the entrance of the reactor, since complete conversion of the reactants is assured. If the hot spot temperature is too hot, though, the catalyst will slowly sinter and deactivate. This dichotomy of goals can be reconciled using automatic control of the hot spot temperature. The control objective is to regulate the hot spot temperature for a variety of different operating conditions including different temperature setpoints and inlet concentrations. Regulation of the hot spot temperature is maintained by varying the N2 flow rate: large flow rates are necessary to keep the hot spot temperature from exceeding a preset upper bound, while lower flow rates are needed to keep the hot spot inside the reactor. Since the product gases, methane and water, act as inerts, control of the nitrogen flow rate is analogous to control of a mass recycle stream around the reactor. In this control configuration, the product gases are recompressed and recycled back to the reactor inlet. Increasing the flow rate through the reactor reduces the
Ind. Eng. Chem. Res., Vol. 34, No. 5, 1995 1751 370i ,0°1
350;. 2
,
1
I
6
.
I
I
I
I
1
10
I
14
I
.
I
I
i
18
N2 Flow Rate (slpm)
320
1 , 0
Figure 7. Estimated hot spot temperature for a small ramp in Nz flow rate.
temperature of the hot spot whereas lowering the flow rate reduces utility costs. This analogy holds because there is complete conversion inside the reactor, the product gases do not affect the reaction kinetics, and the temperature of the exit stream is approximately equal to the inlet temperature. In summary, the problem under study is to control the hot spot temperature by manipulating the nitrogen flow rate in the presence of disturbances in inlet reactant concentration. The axial temperature profile along the reactor is measured using 16 fured thermocouples along the central axis. These thermocouples are spaced at irregular intervals and are sampled once every 5 s. The hot spot temperature is estimated from these thermocouples by fitting a parabola through the three largest temperatures. The maximum value for the parabola is then computed and used as the hot spot temperature estimate. Likewise, the hot spot position is estimated from the position of the parabolic maximum. Figure 7 shows the response of the estimated hot spot temperature t o a slow ramp in Nz flow rate. While it is believed that the true hot spot temperature is a monotonic function of Nz flow rate, the estimated hot spot exhibits several local minima. These minima, which are a consequence of the hot spot interpolation, place severe limitations on the robust controller because the steady state gain of the local linearized model changes sign as the hot spot passes through a minimum. To avoid these problems, the transfer function relating the hot spot temperature to the Nz flow rate is identified only over the range 10-16 slpm (standard liters per minute) where the interpolation effects are relatively small. While other means of estimating the hot spot temperatures have been considered, parabolic fitting is adopted because of its simplicity and relative accuracy.
7. Application of Relay Controller Using the nitrogen flow rate range specified in the previous section, a relay controller is set up which regulates the hot spot temperature at 335 "C.For this operating range, the crossover frequency of the plant is found t o be 0.016 rads, and the breakdown frequency is approximately 0.17 rads. The hot spot temperature response for these two frequencies is shown in Figure 8. No attempt is made to account for the decrease in the magnitude of the output response which naturally occurs due to the rolloff of the transfer function by increasing the magnitude of the input swing. It is fixed for all experiments.
200
400
600
800
Time (sec) Case A closed-loop response at crossover frequency Case E . closed-loop response near breakdown frequency Case C - open-loop response to a fast square wave input ~
Figure 8. Hot spot temperature response for reactor under closedloop relay control. Table 3. Flow Rate Ranges Used in the Reactor Identification Experiments flow rate steady state hot spot temperature case (slpm) for mean flow rate ("C) A 10-12 362 B 12-14 338 C 14-16 327
To determine whether the breakdown is due to noise or uncertainty, the control loop is opened, and a square wave of frequency 0.28 r a d s is directly introduced into the reactor. The result of this experiment, Figure 8, shows the hot spot temperature to be periodic even though its closed-loop counterpart is not. This implies that the breakdown of the relay experiment is due to noise and not plant uncertainty. This is not surprising since the amplitude of the output is so small.
8. Open-Loop Identification Using a frequency range with an upper limit specified by the breakdown frequency, a PRBS is constructed and implemented for three different flow rate ranges. (See Table 3.) Each binary sequence contains 1023 points and lasts for 5 h. A local transfer function is estimated for each of the three cases in Table 3 using spectral analysis. A global model is then obtained by applying the regions-mapping technique to these local estimates. Figure 9 shows the resulting nonparametric model. As expected, the estimated crossover frequency of 0.014 r a d s is very close to that identified using the relay controller. Surrounding this model are the upper and lower uncertainty bounds. The uncertainty associated with this model is small at low frequencies and remains constant over the range 0.001-0.02 radls. At 0.02 radl s, the uncertainty increases in size and remains relatively large for higher frequencies. Because of the inherent trade off between robustness and performance, it is expected that the closed-loop bandwidth will be smaller than 0.02 rads. 9. Controller Design and Uncertainty Validation
To check the validity of the uncertainty bound, internal model control, IMC, theory is used to design a
1762 Ind. Eng. Chem. Res., Vol. 34, No. 5, 1995
f
0.001 0.001
1
I
0.01
0.1
Frequency (rad/sec)
200,
1
-200 0.001
1
I 0.01
0.1 0.001
0.1
Frequency (rad/sec)
I
0.01
0.1
Frequency (radlsec)
Figure 9. Identified frequency response of methanation reactor
Figure 10. Robust stability, q&, as a function of frequency for
with uncertainty bounds.
different robustness filters.
robust controller and that controller is then implemented on the reactor. (See Morari and Zafiriou (1989) for details on the design method.) The first step of the design procedure is to approximate the nonparametric model by a parametric model. Only stable and causal parametric models were considered in this approximation step. A three pole-three zero model is chosen to represent the nonparametric frequency response and the parameters of this model are fit using linear least squares. This model is (sampling interval is T,= 5 s) P(z) =
(z (Z
+ 0.681)(~- 1.136)(~- 2.241) - 0.975)(~- O.809)(~+ 0.339)
(42)
Alternatively, an H, optimal technique for model fitting is available (Rivera and Morari, 1987; Helmicki et al., 1989). Using a first-order robustness filter,
(43) and the parametric model, Pk),the controller, C(z,II)is parametrized as
p -1
C(z,II)=
M
f
1 - PAf
=
N(zA D1(Z,A)- D,(z,II)
(44)
where
- 0.975)(z - 0.809)(z +
N(z,A) = 112(1 -
0.339) (45)
D1(z,II)= ( Z
+ 0.681)[(~- 0.880)(~- 0.446)(~(46)
D,(z,II) = 0.393(z
+ 0.681)[(1 - e-TJ')z(z 1.136)(~ - 2.24113 (47)
Because there is a mismatch between the parametric model, P,and the nonparametric model, &, the additive uncertainty, 4, is recalculated using
4 ( w k )= max 2
lc(ok) - P(e'""S)I
(48)
-
0 4002
Figure 11. Closed-loop response of system to step in hot spot temperature setpoint. Filter time constant 1 = 47 s.
from which the multiplicative uncertainty, computed.
C,can be (49)
Figure 10 shows a plot of the complementary sensitivity, 17, times the multiplicative uncertainty for three different 1. For a controller to robustly stabilize the plant, q& must be less than 1 at all frequencies (Morari and Zafiriou, 1989). This is true for all II I 164 s. Therefore, if the calculated uncertainty bound is a tight bound, II = 164 s should stabilize the reactor over the whole operating window, whereas II < 164 s will be closed-loopunstable for at least one operating condition. A number of validation experiments are conducted in which the hot spot temperature setpoint is given either a positive or a negative step change. Several filter parameters in the range 7.2 s I1 I497 s are tested. Initial experiments reveal that it is more difficult to stabilize the system around hotter hot spot temperatures. Therefore, most of the validation experiments are conducted with a final setpoint temperature near 360 "C. It is found that for 1 I47 s the closed-loop system is unstable while for 1 2 164 s the closed-loopsystem is stable. For the intermediate range 47 s < II < 164 s responses resembling limit cycles are obtained and no clear conclusion about stability can be drawn. Typical experimental runs for 1 = 47 s, 1 = 97 s, and 1 = 164 s are shown in Figures 11,12, and 13, respectively. These experimental results tend to validate the uncertainty bounds obtained with the regions-mapping method. In addition to the setpoint changes, the system is also examined for disturbances in the COn flow rate. This
Ind. Eng. Chem. Res., Vol. 34,No. 5, 1995 1753
Time (sec)
3
Figure 12. Closed-loop response of system to step in hot spot temperature setpoint. Filter time constant A = 97 s.
0
1000
O
2000
0
1500
2000 A
Figure 15. Closed-loophot spot response for 15%increase in COz flow rate. Setpoint temperature is 335 "C.
- ,,,]
B 330
1000 Time (sec)
0
I
5003
A
3000
Time (sec)
Figure 13. Closed-loop response of system to step in hot spot temperature setpoint. Filter time constant 1 = 164 s. Time (sec)
331 7
Figure 16. Closed-loop hot spot response for 15% increase in C02 flow rate. Setpoint temperature is 365 "C.
0
500
1000 Time(sec)
1500
2000
Figure 14. Closed-loop hot spot response for +15% increase in COz flow rate. Setpoint temperature is 325 "C.
is an important and often analyzed regulatory problem since the inlet feed concentrations of an industrial reactor are often fmed by upstream conditions over which the control system has no effect. In this reactor, increases in the reactant concentrations are known to increase the hot spot temperature. The goal of the control system is t o maintain the hot spot temperature at its setpoint value in light of these changes. For this set of experiments, the system is initially brought to steady state a t a predefined setpoint. A +E%change in the inlet COz flow rate is made and the closed-loop response of the hot spot temperature is monitored. If left uncorrected, a 15% increase in the C 0 2 flow rate raises the steady state hot spot temperature from 10 to 30 "C depending on the total flow rate of gas through the system. Using the previously designed controller with a filter time constant of 164 s, the closed-loop responses of the hot spot are examined around three different operating points: hot spot temperatures of 325,335,and 365 "C. The results of these experiments are shown in Figures 14-16. From these figures, it is apparent that the controller maintains its robust stability even though it has not been explicitly designed around the new operating conditions. 10. Conclusions
The objective of this work is to design a robust model based controller for an experimental packed bed reactor system. The key difficulty for designing such a control-
ler is to identify reliable uncertainty bounds from experimental data. To decrease the potential conservativeness of this design, nonparametric models and their corresponding bounds are considered. Three uncertainty bounds for this model are derived: two bounds based on the residual spectrum and a third new bound based on a regions-mapping technique. These uncertainty bounds are examined using Fourier analysis, and it is shown that the uncertainty bounds are composed of two parts: one which measures the deviation of the model from the empirical transfer function estimate and a second which measures the noise level. Next, it is shown how a relay controller preceded by a lead compensator can be used to determine the frequency range over which the control identification experiment should be conducted. Relay experiments were previously proposed in the literature as a technique to identify uncertainty bounds. However, for the reactor system, it was found that this relay experiment does not help in locating the uncertainty bounds since the high-frequency response is dominated by noise. Then a simple nonlinear model is used to test the accuracy of the various identification methods. All three uncertainty bounds are compared to theoretical bounds obtained for the linearized model, and the bounds based on the residual spectrum are discarded as being deficient. This result is especially important for systems which need to be identified over a wide range of operating conditions. Finally, the identification of a fixed bed reactor is performed and the uncertainty bound using the regionsmapping technique is estimated. The transfer function sought is the change in hot spot temperature to a change in inert gas, N2, flow rate. Using the identified model and the computed uncertainty bound, an IMC controller is designed which is shown t o be robustly stable if the time constant of the robustness filter, A, is greater than 164 s. This agreed with experimental results which showed the closed-loop system to be stable for 1 2 164 s, but unstable for A I47 s. Although the uncertainty description is based on a fixed reactant flow rate, the
1764 Ind. Eng. Chem. Res., Vol. 34,No. 5,1995
robust controller effectively regulates the hot spot temperature for a +E%change in the reactant flow rate.
Acknowledgment H.M.B. acknowledges the support of the Rothschild and Bantrell Foundations. Literature Cited Abramowitz, M.; Stegun, I. Handbook of Mathematical Functions; Dover Publications: New York, 1964. Astrom, K.; Hagglund, T. Automatic tuning of simple regulators. Proceeding from ZFAC Ninth World Congress, Budapest; Pergamon Press: Oxford, 1984; p 1867. Bartlett, M. Smoothing periodograms from time-series with continuous spectra. Nature (London) 1948,161,686. Blackman, R.; Tukey, J. The Measurement of Power Spectra; Dover Publications: New York, 1959. Brillinger, D. Time Series-Data Analysis and Theory; Holden Day: San Francisco, CA, 1975. Cutler, C.; Ramaker, B. Dynamic matrix control-a computer control algorithm. Proceedings Joint Automatic Control Conference, San Francisco, CA American Automatic Control Council: 1980; Paper -5-B. Daniell, P. Discussion following “On the theoretical specification and sampling properties of autocorrelated time series” by M. S. Bartlett. J. R. Stat. SOC. 1946,8, 88. Gardner, W. Statistical Spectral Analysis-A Nonprobabilistic Theory; Prentice-Hall: Englewood Cliffs, NJ, 1988. Godfrey, K. Theory of the correlation method of dynamics analysis and its applications. Meas. Control 1969,2,65. Goodman, N. Statistical analysis based on a certain multivariate complex Gaussian distribution (an introduction). Ann. Math. Stat. 1983,34,152. Hannan, E. Multiple Time Series; John Wiley and Sons: New York, 1970. Helmicki, A,; Jacobson, C.; Nett, C. H, identification of stable lsi systems: A scheme with direct application to controller design. Proceedings from the American Control Conference, Pittsburgh, P A Department of EECS, Northwestern University: Evanston, IL, 1989; p 1428. Jenkins, G.; Watts, D. Spectral Analysis and Its Applications; Holden-Day: San Francisco, CA, 1968. Kay, S.; Marple, S. Spectrum analysis-a modern perspective. Proc. IEEE 1981,69,1380.
Koopmans, L. The Spectral Analysis of Time Series; Academic Press: London, 1974. Laughlin, D.; Morari, M. Graphical stability analysis for control systems with model parameter uncertainties. Comput. Vision Graphics Image Process. 1989,47,59. 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. Ljung, L. On the estimation of transfer functions. Automatica 1986,21,677. Ljung, L. System Identification-Theory for the User; PrenticeHall: Englewood Cliffs, NJ, 1987. Loh, A.; Corrba, G.; Postlethwaite, I. Estimation of uncertainty bounds for robustness analysis. ZEE Proc., Part D 1987,134, 9. Morari, M.; Zafiriou, E. Robust Process Control; Prentice Hall: Englewood Cliffs, NJ, 1989. Palazoglu, A.; Owens, S. Robustness analysis of a ked-bed tubular reactor: Impact of modeling decisions. Chem. Eng. Commun. 1987,59,213. Postlethwaite, I.; OYoung, S.; Gu, D.; Hope, J. H” control system design: A critical assessment based on industrial applications. In Tenth ZFAC Triennial World Congress, Munich Pergamon Press: London, 1987; p 301. Rivera, D.; Morari, M. Control-relevant model reduction problems and p-controller synthesis. Znt. J. Control for SISO H2,L, 1987,46,505. Rivera, D.; Webb, C.; Morari, M. A control-relevant identification methodology. In Proceedings from AIChE Annual Meeting, New York; AIChE: New York, 1987. Smith, R. Personal communications, 1989. Smith, R.; Doyle, J . Model invalidation: A connection between robust control and identification. Proceedings from American Control Conference, Pittsburgh, P A Department of EECS, Northwestern University: Evanston, IL, 1989; p 1435. Vidyasagar, M. Nonlinear Systems Analysis; Prentice-Hall: Englewood Cliffs, NJ, 1978. Webb, C. Robust Control Strategies for a Fixed Bed Chemical Reactor. Ph.D. Thesis, California Institute of Technology, 1990. Received for review July 5, 1994 Revised manuscript received February 2, 1995 Accepted February 16,1995@ IE940409V Abstract published in Advance ACS Abstracts, April 1, 1995. @
