2656
I n d . Eng. Chem. Res. 1994,33,2656-2667
Experimental Design Considerations for Dynamic Systems Roger W. Shirt,?Thomas J. Harris,’ and David W. Bacon Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada K7L 3N6
Although experimental designs have been extensively applied for steady-state modeling of process behavior, their use in studying dynamic systems has been more limited. Optimal designs were primarily focused on minimizing uncertainty in the transfer function parameter estimates, until Ljung pointed out the importance of minimizing bias in the assumed form of the transfer function. This paper reviews existing methodologies for designing input test signals and draws comparisons among them through the use of simulation examples. Results indicate t h a t input designs which emphasize relevant frequencies for the controller application give superior performance when either an adequate process transfer function model has been specified or bias in its form exists. In the latter case, minimization of the frequency domain bias is found to be the most important factor in ensuring reliable controller performance. Distributional results for the closed-loop performance criterion provide important insights into the dominant source of model uncertainty.
Introduction Modern model-based control schemes require reliable process models to achieve optimal performance. Often, some form of perturbation signal is applied to a system to provide information for model building or evaluation. Selection of a perturbation signal for identification of process dynamics has traditionally been determined by the particular parameter estimation method used, or alternatively has been made on an ad hoc basis. Common choices of input signals used are sinusoids, impulse and step changes, and pseudorandom binary sequences (PRBS’s). The benefits of statistically designed experiments for characterizing steady-state process behavior have been recognized in the chemical process industry for at least 40 years (e.g., Box and Draper, 1987; Rippin, 1988). Typical objectives of such experiments include (a) detecting the more promising model forms from among a proposed group of candidate forms (the model discrimination problem), (b) providing improved precision of the parameter estimates in one specified model form (the parameter precision problem), (c) guarding against possible inadequacies in a proposed model form (the model bias problem), and (d) confirming process performance as predicted by a specified model, including comparisonswith predictions from alternative plausible models (the model verification problem). Developments in experimental designs for steadystate models have included both linear and nonlinear model forms, and applications have been reported in all fields of engineering and science, as well as in the social sciences. Extension of experimental design concepts to dynamic models has been slow t o develop for several reasons, including (a) the broad variety of forms of process dynamics that can occur, (b) limited familiarity among dynamic modelers with the extensive statistical literature that exists for experimental designs in the steadystate case, (c) difficulty in characterizing uncertainty when an inadequate model form is employed, and (d) the broad variety of possible end uses of a dynamic model for process control. ~
* Author to
~
~~
~~~~
whom correspondence should be addressed.
E-mail:
[email protected]. Current address: Pulp and Paper Centre, The University of British Columbia, 2385 East Mall, Vancouver, BC, Canada V6T 124.
’
However, the potential benefits of experimental design for dynamic systems have been recognized for some time. In 1977, for example, Goodwin and Payne observed that “Until recently most interest in identification was directed towards the data analysis problem. However, there now exists a growing recognition of the importance of experimental design.” Published work during the past two decades has included optimal input signal selection to identify dynamic models in both differential equation forms (e.g., Espie and Macchietto, 1989)and linear discrete time models (e.g.,Mehra, 1974; Goodwin and Payne, 1977; Zarrop, 1979; Ljung, 1987). Prior to the work of Ljung and co-workers during the 1980s, efforts had been directed almost exclusively to the parameter precision problem and the designs employed were usually based upon one or more of the “alphabetic optimal” criteria proposed by Kiefer and Wolfowitz (1959), which are described in the next section of this paper. Each of these criteria involves minimization of some scalar measure of the parameter variancecovariance matrix. It is assumed that an adequate model structure and reasonably accurate parameter estimates are known a priori. In a significant departure from this overall approach, Ljung (1987) proposed an experimental design criterion for linear discrete time models which involves minimization of a squared measure of the uncertainty in the estimated transfer function, accounting for uncertainties in both the selected form of transfer function and the estimates of the parameters in that function. This approach has been extended to continuous time systems by Goodwin et al. (1991). There is now increasing interest in characterizing the transmission of transfer function uncertainty to the resulting control loop performance (Gevers, 1991; Rivera et al., 1990). Results published t o date have normally been limited to linear dynamic models with additive stochastic disturbances. While numerous heuristic and optimal input design methodologies have been proposed, very few direct comparisons have been made among them. Exceptions include Gupta et al. (1976), Zarrop (19791, and Soderstrom and Stoica (1983), who have presented results comparing different alphabetic optimal designs for the case of no model mismatch. The main purpose of this paper is to compare several different input signal designs in their ability to provide good controller performance. PRBS, alphabetic optimal and transfer
0888-5885/94/2633-2656$04.50/0 0 1994 American Chemical Society
Ind. Eng. Chem. Res., Vol. 33, No. 11, 1994 2657 function based designs, are used in different modeling scenarios. The following sections of this paper summarize existing methods for designing input signals for dynamic systems, with special emphasis on the contributions of Ljung (1987). The focus of this summary is on linear discrete time models which are commonly used in modern process control applications. Various design choices are highlighted, reflecting the objective(s)of the experiment. Limitations of these methodologies are noted. A frequency domain approach is shown to be useful in providing insights into the choice of an appropriate experimental design. A simple simulation example is used to compare the effectiveness of five different designs with respect to several performance criteria.
Techniques for calculating M are discussed by Goodwin and Payne (1977) and Zarrop (1979). In the case of a purely stochastic process see Box and Jenkins (1976) or Whittle (19511, who introduced the frequency response equivalent of M. An optimal input signal is chosen to minimize some measure of M-l. Typical measures include the following “alphabetic” criteria (Kiefer and Wolfowitz, 1959):
Experimental Design Criteria
W is a weighting matrix reflecting the relative importance of each component of M-l (9)
A single inputlsingle output (SISO) linear time invariant dynamic system can be represented by a model of the form
where Yt is the measured system output at time t , Ut-l is the (manipulated) input signal at time t - 1,G&) is the process transfer function, z is the unit forward shift operator (i.e., zXt = Xt+l), and Dt is an additive system disturbance at time t. The disturbance Dt is assumed to be a realization of a stationary stochastic process driven by white noise:
D,= H,(z)e,
(2)
where H&) is the disturbance transfer function and et is a white noise sequence. Ho(z)is stable and invertible; i.e., l/Ho(z) is also stable. Let the vector of the P parameters associated with Go(z)and H&) be denoted by 8 and denote estimates of 8 obtained from Npbservations of input-output data from the process by 8. Characteristics of the precision of these parameter estimates gre contained inAthe variance-covariance matrix Cov[8]. Because Cov[OI for dynamic models cannot, in general, be evaluated before input-outpu! data have been collected, an approximation for Cov[81 is required. The Cramer-Rao inequality states that for unbiased estimators of the parameters (Whittle, 1951; Goodwin and Payne, 19771, Cov[hl L M-’
.(3)
where the inequality is understood to imply a matrix norm and M denotes Fisher’s Information Matrix, defined by
where EyIe is the statistical expectation of Y conditional on 8 and Z(Yl8) is the likelihood function. Equality is approached in expression 3 as N -. Thus, for a sufficiently Large sample size, M-’ is a good approximation of Cov[BI. For large sample sizes the average per sample information matrix M is usually used, where
-
M = lim &)M N-
(5)
Det[M-l], a D-optimal design
(6)
Tr[W-M-l], a n A-optimal design
(7)
A,,[M-’I,
(8)
a n E-optimal design
and Det[ I, Tr[ I, and[,A I are the determinant, trace, and maximum eigenvalue, respectively, of the matrix [ I. Minimization is carried out subject to constraints on the variances andor amplitudes of the input and the output. The resulting input signal is usually formulated as a sum of sinusoids (Goodwin and Payne, 1977). Recently, Tulleken (1990) proposed a method for constructing binary test signals which are D-optimal. Because of the nonlinearity (in the statistical sense) of the dynamic model forms, the elements of M-l depend upon the true values of the parameters 8, which are unknown. The “optimality” of a design therefore depends upon the quality of the a priori estimates of 8, as well as on other characteristics (Box and Draper, 1987). For dynamic system models it is important to note that optimal input designs cannot directly influence the precision of the parameters which are uniquely associated with the disturbance component of the model. However, as Zarrop (1979)has shown, those parameters do influence optimal designs for the parameters associated with the process transfer function. As stated above, the objective of an alphabetic optimal design is to minimize some measure of parameter uncertainty. This objective is not specifically related to the end use of the model unless improved parameter precision happens to be the ultimate objective. Watts and Minich (1972) used alphabetic optimal designs in a time domain setting to obtain input signals for minimum variance control. In the next section of this paper it will be seen how the transfer function based input design procedure proposed by Ljung (1987) leads to a logical choice of frequency weighting for variance optimal designs, reflecting the model end use.
Transfer Function Based Input Signal Design Employing a frequency domain approach, Ljung (1987) has proposed a comprehensive methodology for designing input signals. His design criterion accounts for deficiency in the proposed form for the system transfer function, uncertainty in the parameters in that transfer function, and the model end use. The alphabetic optimal design criteria described in the previous section are subsets of this more general criterion. A frequency domain formulation allows further insight into the choice of an appropriate optimal design. For example, the frequency distributions for both model bias and parameter uncertainty can be explored, as well as the effects of various factors (e.g., prefilters, input
2658 Ind. Eng. Chem. Res., Vol. 33, No. 11, 1994
signal) on these distributions. Ljung (1987) and Goodwin et al. (1991) provide further discussions of these features. Some definitions are required to provide a framework for Ljung's criterion. The true system (i.e., the actual dynamic behavior of the process) is described by eqs 1 and 2. Note that the true system is assumed to be linear in the inputs, an assumption which is only approximately true in practice. The fitted process model is
+
Y , = B(z,e)v,-,
H(Z,&e,
(10)
where G(z,& and fi(z,& are the fitted model forms representing the transfer function and disturbance model components, respectively, using the (not yet obtained) parameter estimates 8. Define the following vectors of transfer functions: To(z) = [Go(z) Ho(z)lT
(11)
T(z,e) = [ G ( Z , P ) H(z,8)IT
(12)
The frequency domain transfer function bias can then be defined as T(ei",b) = T(eio,e)- TO(e'")
(13)
where eio has been substituted for z , and --n < w < -n. An input design cost criterion can be defined as
where C(w)is a 2 x 2 user specified frequency weighting function and E&is the expected value with respect to 8. J can be interpreted as a weighted frequency domain quadratic measure of the bias between the true and estimated system transfer functions. Unlike the alphabetic optimal input signal designs, there is no restriction that T(z,8) must have the same structural representation as To(z). As with alphabetic optimal designs, minimization of the design cost function may be carried out subject to constraints on the input and output signals. The optimal input signal can be characterized in terms of its spectrum, @$Pt(w). The input design problem is then minimize J wrt @u
subject to G
(15)
where Gis the set of all applicable constraints. Ljung (1987) has shown that, asymptotically, J can be split into two distinct components, representing contributions from imprecision in the parameter estimates (variance) and inadequacy in the assumed forms of the transfer functions (bias). That is,
J XJ v +
J B
(16)
The individual subproblems of eq 15, minimize J wrt @u
B
minimize Jv wrt @u
subject to G
(17)
subject to G
(18)
are of interest when either JB or JV dominates.
The relative importance of different frequencies is dictated by the end use of the model and is reflected in the input design criterion 14 by the frequency weighting function C(w). That is, the variance and bias are assessed with respect to the intended end use of the fitted model. Steady-state designs also employ a weighting function over the operating range of interest (Box and Draper, 1987). Appropriate forms of C(o)have been derived by Ljung (1987)for process simulation, one-stepahead prediction, minimum variance control, and pole placement regulation. C(o) usually contains a frequency description of the expected input Ut once the model has been implemented (possibly in a closed-loop control scheme). Unless a completely new control application is being considered, this component of C(w) can be estimated from existing plant operating data. Bias Optimal Designs. Let 8* denote the asymptotic estimates of the parameters 8. Then 8* can be regarded as the best estimates of the parameters for the chosen model foym. One definition of transfer function bias is then [G(eiw,8*)- Go(eiw)l,the frequency domain difference between the best estimate of the assumed transfer function model form and the true system transfer function. Clearly, if the selected model form is (or includes) the true model form, then the bias will be zero. If a reliable estimate of the process transfer function alone is of interest, then the bias contribution t o the total design cost, J B , is (Ljung, 1987)
Go(e-iw)lCll(w)dw (19) where C l l ( w ) is the (1,l)element of C ( w ) . For this form of J B , the solution to problem 17, subject to an input power constraint, is
where a is a constant chosen to satisfy the input power constraint. The superscript on 8 has been suppressed to indicate that the disturbance model form and the parameter values in that form have been fixed a priori for the input signal design stage. Although the disturbance transfer function does not appear explicitly in eq 19, it does affect 8* and in turn the optimal input signal given by eq 20. Knowledge of the disturbance spectrum is always required for calculation of an optimal input signal. Variance Optimal Designs. If an estimate of only the process transfer function is of interest, the expression for JVis given (Ljung, 1987) by
where G is the (1 x P ) matrix of derivatives of G with respect to 8 evaluated at 8*. Thus, minimizing JV is equivalent to minimizing a weighted trace of the parameter variance-covariance matrix. If Cov[BlN is approximated by the inverse of Fisher's Information Matrix, this criterion is equivalent to that for Aoptimality. The weighting matrix W in eq 7 is given by the integral portion of eq 21. The variance cost function (eq 21) is then seen to be a special case of the overall cost function (eq 14) in which the A-optimal
Ind. Eng. Chem. Res., Vol. 33, No. 11,1994 2669 design criterion is used under the assumption that there is no uncertainty in the form of the estimating model. However, in both steady-state and dynamic modeling applications, D-optimal experimental designs are the dominant choice for minimizing variance because of the relative ease with which parameter constraints can be accommodated and the availability of software packages for implementing this design criterion. Also, the information matrix does not need t o be inverted as is the case for A-optimal designs. Ljung has also proposed an alternative method for solving problem 18 based on the assumption that if the model order of the estimating process transfer function is suficiently large, it will provide an adequate representation of G(eio,6) for input design purposes. Although estimation problems can arise from an overparametrized model, this approach leads to an input design solution that is more tractable than those using alphabetic optimal designs, especially for higher order systems. For an input power constrained design obtained under the above assumption, the optimal input spectrum is given by
(22) where @d(O) is the spectrum of the disturbance and ,u is adjusted to satisfy the power constraint. Again we see that even though the focus of interest is the process transfer function, the optimal design depends largely on the spectrum of the disturbance. Ljung has also developed a solution to the more general problem of minimizing the overall cost criterion given by eq 14. Once again an infinite order approximation for the process transfer function is employed and the solution is equivalent to eq 22 above. The equivalence of these two solutions is not surprising. As the process model order tends to infinity, the bias in the transfer function estimate is removed, transforming the problem into one of minimizing the variance contribution.
Implementation of Input Design Methods To illustrate some of these design methods, let us consider the case for which the true process is a second order system with additive disturbance:
where 61 = 1.7236, 82 = -0.7408, and the gaing = 1. The step response of this process is shown in Figure 1. The true stochastic disturbance Dt follows a first order autoregressive model
D,= 0.9Dt-,
+ e,
(24)
where {et} is a normal white noise stochastic sequence with mean value zero. The variance of Dt was fixed a t 1.0 by specifying the variance of {et} to be 0.19. Reasons for using this particular model as an example are the following: (i) the parameters are close to the stability boundaries for a discrete second order system, (ii)there is little response in the first sampling interval to changes in the manipulated variable, and (iii) the
1
.
5
~
1
0
0
10
30
20
40
50
Time Figure 1. Step response.
stochastic model represents drifting, meandering disturbances which are frequently encountered in practice. Open-loop identification experiments for the above system were simulated. Five experimental designs (given later) for the process input were evaluated with respect to the resulting parameter uncertainty and the impact of using the estimated model in a feedback control scheme. The parameter uncertainty was measured by (i)the trace and determinant of the variancecovariance matrix of the estimated system parameters and (ii) the joint confidence region of the real and imaginary components of the Nyquist plot at three frequencies. The impact of using parameter estimates for controller design depends on the controller design method. A model following controller can be designed so that the closed-loop dynamics are the same as the open-loop dynamics of the estimated model (Astrom and Wittenmark, 1984). Define the performance degradation, ft, to be the difference at time t between controller responses based on the estimated and true process models. Then N
APerf =
eft" t=l
is a measure of model performance in the control scheme. We note that this measure of closed-loop performance is equivalent to that given by Ljung (1987). We examined APerf for a unit setpoint change to the model following controller designed using the estimated transfer function. The designs were evaluated against the above criteria in three scenarios: (i) model with no bias, (ii) undermodeling (plant transfer function represented by a first order model with deadtime), and (iii) using extended prediction horizons for parameter estimation. The input designs were the following: M1, Pseudorandom binary sequence (PRBS) with a switching time of 5 sampling intervals; M2, PRBS with a switching interval of 10 sampling intervals; M3, bias optimal design for control (eq 20); M4, variance optimal design for control (eq 22); M5, D-optimal design for parameter estimation (eq 6). The switching times for designs M1 and M2 were randomly selected. For designs M3, M4, and M5 the spectrum of the input signal was determined using the true process transfer function for the zero bias case, and a priori values based on an approximate first order fit in the case of a biased model. Both scenarios used the true disturbance transfer function. To obtain a time domain realization of the M3 and M4 inputs in this study, a low order autoregressive
-
2660 Ind. Eng. Chem. Res., Vol. 33, No. 11,1994
moving average time series model of the following form was used first as an approximation:
-,+...+
Ut=AIUf_,+...+ ApUtI,,+at+Bla,
Bqa,-, (26) where ofrepresents the departure of Ut from its mean value and {a,} is a white noise stochastic process with zero mean. The parameters {A;, Bi} in eq 26 were adjusted t o provide an input series whose spectrum closely matched the desired input spectrum. Next,_for input signals M3, M4, and M5, the input sequence {UT} was transformed into a binary sequence by clipping the individual values (Kedem, 1980; Ljung, 1987):
This binary sequence was then used as the input signal. The process defined by eqs 23-27 was simulated for
t = 1,2, ..., 500. The parameters of the transfer function
(first or second order model) for the system were estimated using a k-step-ahead prediction error criterion with a constrained optimizer to ensure that the parameter estimates remained within the stability region and that the process gain was within the arbitrarily selected interval 10.5, 2.51. This approach presumes accurate prior knowledge of the process, and its advantages are discussed at length by Tulleken (1993). Tulleken implemented the constrained optimizer by imposing n simultaneous inequality constraints on the 8 s . where n is the order of the system. When identifying a second order model we chose to reparametrize b(z) as follows:
61
y1=1-6,
(28)
Y z = 6,
(29)
C""(S) Nyquist EIIipm
Figure 2. Simulation procedure
With this representation, stability of the polynomial b(z) is ensured when the y's are each bounded by f l . This parametrization is equivalent to using the partial correlation representation of an autoregressive relationship (Box and Jenkins, 1976). Extensions to higher order systems are straightforward. The real and imaginary components of the estimated open-loop transfer function were calculated at those frequencies where the amplitude ratio of the true process was 0.9, 0.707, and 0.1. The model following controller was then implemented and APerf for a unit set point change was calculated. This entire procedure (input realization parameter estimation controller simulation) was repeated 1000 times for inputs amplitudes 2.0, 2.5, and 3.0 and is summarized in Figure 2. The results of the simulations were summarized as follows. First, the iterative procedure outlined by Rousseeuw and Leroy (1987) was used to calculate the parameter variance-covariance matrix. From this, approximate 95% confidence ellipses for the true real and imaginary components of Go(e'O) at frequencies corresponding to amplitude ratios of 0.9,0.707, and 0.1 were plotted on the Nyquist diagram. In the case of no model mismatch, the required assumption that these components follow a bivariate normal distribution can be justified using asymptotic theory as follows. Denoting the real and imaginary components of the gradient of the system transfer function with respect to the
+
+
system parameters as (VOGO), and (VeGoX, respectively, the real and imaginary components of the system transfer function can be expanded in first order Taylor series about their nominal values as
N,= (VOGo)+O Ni = (V,G,),AO
(30)
If e(t) in eq 2 is normally distributed, then so are the parameter estimates with variance-covariance matrix M-l. Asymptotically then, the real and imaginary components also follow a bivariate normal distribution, centered about their nominal values with variancecovariance matrix:
Recently several authors have presented results for producing elliptical Nyquist regions which encompass both parameter and model uncertainty (Gevers, 1991; Wahlberg and Ljung, 1992; Goodwin et al., 1992). In this paper we plot only the variance components as our interest is in making qualitative statements about the relative magnitudes of frequency domain bias and variance. Here the size of an ellipse is related to the
Ind. Eng. Chem. Res., Vol. 33, No. 11,1994 2661 uncertainty in the real and imaginary components and its position in relation to the true Nyquist curve indicates the degree of bias a t that particular frequency. An asymptotic expression can also be developed for the variation in closed-loop performance as follows. Assuming a fimt order expansion about the true process model to be adequate, the derivative of the closed-loop tracking error is
0
Frequency (rad) where ysp denotes the set point value. Evaluate this expression a t 8 = 8*, the converged parameter values. The controller transfer function, G,, is based on the estimated model. When the model bias is negligible, it can be shown that an approximation for the variation in performance degradation is
-
= 1/(2n) AT8fFT(e-J") F(e'")do A0
APerf= t=o
(33)
where A 8 = 8 - 80, which is asymptotically normally distributed. Consequently, the variable APerf is distributed as a weighted sum of chi-squared variables (Box, 1954) with
0
li
Frequency (rad) Figure 3. (a) Theoretical input spectrum. (b) Realized input spectrum
mean{APerf) = trace(WM-l) VarIAPerf) = 2{trace(WM-')'}
(34)
where W denotes the value of the contour integral. Exact percentage pointe and approximations for weighted sums of chi-squared variables are discussed in MacGregor and H a m s (1993) and H a m s et al. (1993) and the references contained therein. The consequence of this result is that input signal designs which produce better average controller performance also produce tighter deviations in this performance when the conditions under which eq 34 was developed hold. In summary, the design variables found to influence the results from these simulations were the input spectrum, the magnitude of the input signal, the model structure, and the estimation method. The response variables were the performance measures for parameter uncertainty and closed-loop control. The Zero Bias Scenario. Correct model structures for the process and disturbance were employed for estimation in this case. The input spectra and resulting time domain representations from designs M3, M4, and M5 are summarized in the Appendix. Figure 3a shows the theoretical spectra of the manipulated input variables for each of the five designs considered. For example, the bias optimal design (eq 20) tends to emphasize lower frequencies than does the variance optimal design (eq 22). Figure 3b shows the averaged unsmoothed input signal periodograms computed over the 1000 realizations window. For the case of the D-optimal input (M5) the desired spectral features are not matched very well. This is mainly due to the fact that the lower of the sinusoidal frequencies required for the input realization has a period of 1851 sampling periods compared with a total experimentation time of 500 sampling periods. As well, note that the M5 spectrum has only one dominant frequency when support a t a t least two is required for persistence of excitation. This may lead to poor condi-
2
3 4 Input Signal Amplitude
5
Figure 4. Trace of normalized covariance matrix, second order model.
tioning of Cov(8). In contrast, the agreement of each of the other input signals to the theoretical spectra is generally good. As binary input signals were used in all simulations, we can partition the input spectrum into amplitude and frequency components Q,(w) = A2&,,(o)
(35)
The frequency component is based on unit amplitude. Figure 4 shows the relationship between input signal amplitude and trace of the variance-covariance matrix for the zero bias scenario. The expected inverse relationship with amplitude is immediately apparent. It is also evident in each case that a family of curves is obtained based on the frequency characteristics of the input signal. Due to the poor realization of the Doptimal input, the parameter precision is relatively low when compared with that from any of the other methods. In general, the input signals which emphasize lower frequencies provide greater parameter precision for this modeling situation.
2662 Ind. Eng. Chem. Res., Vol. 33, No. 11, 1994
O 2
I
M3
u=2 5 I
1
{
0.1
2
3 4 Input Signal Amplitude
5
2
3 4 Input Signal Amplitude
5
u=2 5 1
0
L
02
-
00
0 2
1
I
04
06
I
08
10
12
Real Component Figure 5. 95% confidence regions, second order model.
Maximizing the signal-to-noise (STN) ratio at important frequencies has been used to motivate input designs (Ljung, 1987; Goodwin et al., 1991; Rivera, 19921, including eq 20. The frequency domain interpretation of the STN is
A frequency domain expression for the information matrix has been given by Ljung (1987):
1.5
1 I
fa
0
l 0
a
",
v
'0,
Each expression is composed of a squared amplitude contribution and a frequency weighted contribution from the input signal. However, the frsquency weightings are different, most notably with M influenced by the shape of the disturbance spectrum. Thus the overall STN should only be used as a rough indicator of input signal quality. Approximate 95% ellipsoids for the real and imaginary components of the process transfer functions resulting from input designs M1 and M3 are shown in Figure 5. Note that the ellipsoids, which are centered a t the true Nyquist frequencies, shrink with increasing input amplitude. Differences in size, which parallel those observed in the variance-covariance matrix results, are also observed due to input signal frequency content. The orientations of the ellipses also change with frequency and input signal frequency content, but not with input signal amplitude. These observations are consistent with the asymptotic theory results of eqs 31 and 35. Figure 6a reveals the relationship between Mean(APer0 and input signal amplitude to be strong for all inputs. Significant differences in controller performance result from the different input signal frequency contents with lower frequency input signals generally producing improved results. These results parallel those of Figure 4 for precise parameter estimation. Std(APerf) results (Figure 6b) are similar t o those for Mean(APer0. A plot of Std(APerf) against Mean(APerf)(omitting M5 values) is shown in Figure 6c, and it reveals a straight line relationship with slope of 1.7. This is in good agreement with the theoretical value of
uMI
(I)
v+
8
OM2
M3 M4 MS
d"
0-
7
0
0.5
1
1.5
Mean(aPeQ Figure 6. (a) Closed-loop performance measure, second order model. (b) Variation in closed-loop performance measure, second order model. ( c ) Closed-loop performance measures, second order model.
70
Mean(aPerf) Figure 7. Distribution of closed-loop performance measure, second order model.
1.4 computed from eq 34. M5 values appear not to follow this Telationship, possibly due to poor conditioning of Cov(6). Figure 7 shows the corresponding distribution of APerf values for M5 to be much flatter than that for M3.
Ind. Eng. Chem. Res., Vol. 33, No. 11, 1994 2663
0.8
.
5
I
1.4
I
!
1
I
1.2
1
0.2 0
4.-
r
2
3
I
1 04 4
5
2
Input Signal Amplitude
,
3
7
-
5
4
Input Signal Amplitude
Figure 8. Trace of normalized covariance matrix, first order model.
1.5
1
The Biased Model Scenario. The (second order) transfer function model was approximated using a first order model with dead time:
G, =
g(1 - d ) f b
1 - 6z-l
(38)
The disturbance model was first order autoregressive. Optimal input designs require a priori specification of the parameters of the nominal model. For this case the values used were g = 1, 6 = 0.95, and b = 3, producing a step response for the approximate model that closely matched that of the true model. We note that the D-optimal input was a single sinusoid of frequency w = 0.068 with a period of 92 sampling intervals. The simulation procedure was identical to that described for the zero bias scenario. The time delay, b, was estimated jointly with the system and disturbance parameters using a one-step-ahead prediction error criterion. The parameters were normalized with respect to their average values for the entire set of simulations. The controller was designed on the assumption that the estimated first order transfer function model was correct. Tr(Cov(8)) is dominated by uncertainty in the dead time parameter, b, and the results in Figure 8 reflect this fact. Results are nearly opposite those from the zero bias scenario, with higher frequency test signals giving better parameter precision. Conversely, Figure 9a shows that controller performance improves with lower frequency input test signals. Overall performance is worse than in the second order model case and larger differences are seen between methods. We note that the bias optimal input provides the best performance among the input design methods considered. The behavior of Std(APer0 does not closely follow that of Mean(APer0, as shown by Figure 9b. Input signal amplitude does not always lead to a reduction in Std(APerf), particularly for higher frequency inputs. Clearly this is an undesirable feature of an input test signal. A n unweighted least squares regression of Std(APerf) against Mean(APer0 gives a slope of 0.7 (Figure 9c), well below the value 1.7 obtained for the case with no model mismatch. These results can be given a frequency domain interpretation by considering Figure loa, which shows Nyquist ellipses for M1, M3, and M4 inputs with amplitudes of 2.5. Generally, the ellipses are not centered over the true Nyquist curve due to estimation bias. When the location of any ellipse converges toward the nominal Nyquist curve, the behavior of Std(APer0 matches more closely that of Mean(APerf). Figure 10b
I 0
1 2
1.5
fa
-
7
3 4 Input Signal Amplitude
5
1
I
1 -
d
I
v M4
M5
I
0 0
0.5
1
1.5
Mean(APerf) Figure 9. (a) Closed-loop performance measure, first order model. (b) Variation in closed-loop performance measure, first order model. (c) Closed-loop performance measures, first order model.
demonstrates that both the size and orientations of the ellipses are affected by input signal amplitude. These figures clearly demonstrate that bias, and not precision, a t any particular frequency governs the behavior of APerf. The Biased Model: Extended Prediction Horizon Scenario. In this scenario the parameter estimates were those values that minimized the sum of squares of the five-step-ahead prediction errors of the first order model given in eq 38. Input designs remained the same as in the one-step-ahead prediction scenario. Comparison of Mean(APerf)for this case (Figure l l a ) with Mean(APerf)for one-step-ahead prediction (Figure 9a) shows marked improvements in the performances of the M4 and M2 inputs. The M1 input now has a stronger relationship with amplitude but has sacrificed some performance at lower input amplitudes. The relationship of the behavior of Std(APerf), shown in Figure l l b , with that of Mean(APerf)is more consistent with eq 34, as shown in Figure l l c . The slope of the regression line is 1.6 with a significant nonzero intercept. Overall improvements in behavior are attribut-
2664 Ind. Eng. Chem. Res., Vol. 33, No. 11, 1994 1.4
0.2
1.2 0.0
f 1
c
C a,
$ 0.8
c
v
5 0.6
0
E -o.2
8 0.4
0
0
0.2
5 - 0 4
.-C 0) m
E
0
v-
2 -06
-
, -3 4 5 Input Signal Amplitude
3 -0
a
2.5
fa 2
-1 0
-02
00
02
04
06
08
10
12
Real Component
a 1.5 v B v , 1 0.5
0 c
0.0-
~-
1 -
,7
2
5
3 4 Input Signal Amplitude
C
Q)
C
0 Q
E 0 0
-0.2 -
25
7 ---I -
21
g1.5
r
LL
a
M3 u=2 0
7
5 - 1
UMl
v)
OM2
OM3
0.5 -0 a -
.1 0
M4
M1
rn MS
u=5 0
+ -02
00
0
0.5
1
1.5
2
2.5
Mean(APerf) I
I
1
I
1
02
04
06
08
10
12
Real Component Figure 10. (a) 95% confidence regions, first order model: input signal magnitude i 2 . 5 . (-1 M1;(- - -1 M3; M4.U = 2.5. (b)
Figure 11. (a) Closed-loop performance measure, first order model: extended prediction horizon. (b) Variation in closed-loop performance measure, first order model: extended prediction horizon. (c) Closed-loop performance measures, first order model: extended prediction horizon.
(.e*)
95% confidence regions, first order model.
able to reductions in frequency domain bias, as illustrated by comparing Nyquist ellipse positions in Figure 12 with those in Figure loa. The ellipses in Figure 12 are more centered over the nominal Nyquist curve. However, performance improvements are generally not evident for the M1 and M5 inputs. We note the significant increase in the size of the M1 Nyquist ellipse with the extended horizon scheme and see that precision has been traded off for bias reduction.
Summary and Conclusions All (‘optimal”methods attempt to shape the input signal frequency properties according to some criterion of goodness of model fit or controller performance. Design procedures depend on the assumed model for process dynamics and the degree of confidence in this model, which gives rise to “bias” and ‘(variance”designs. Frequency weighting of the input design may be em-
ployed which reflects the anticipated use of the model (i.e., parameter estimation, controller application). Several measures of input design performance have been employed in this study. In the idealized case of correct choices for the process and disturbance model forms, simulation results indicate that parameter precision is an adequate indicator of controller performance. However, when model mismatch is present, parameter precision is not a good indicator of controller performance and should not be used as such. Nyquist uncertainty regions constructed according to the procedures described under Implementation of Input Design Methods allow one to independently compare the resulting bias and variance a t a particular frequency. It was seen that they have orientations which are frequency dependent and that there is a decrease in size with increasing input signal amplitude and properly selected frequency content. These observations apply for the cases of both correctly and incorrectly specified transfer functions when the bias at that particular
Ind. Eng. Chem. Res., Vol. 33, No. 11, 1994 2665 Table 1. Percentage Improvement of MeadAPerf) and Std(APerf)(in Parentheses) Using M3 (Bias Optimal) Input for Test Signal Amplitude of 2.5 ~~~
modeling scenario zero bias bias k = 1 bias k = 5
0
2 -o.2
sE E
O’Ol
-0.6
-08
-1 0
-02
00
02
04
I
I
06
08
!A -
10
12
Real Component Figure 12. 95% confidence regions, first order model: extended prediction horizon. Input signal magnitude f 2 . 5 . (-) M1; (- - -) M3; (. M4.
-
a)
frequency is minimal. Orientations of the ellipses are affected by input signal amplitude when frequency bias is present. An expression has been developed which shows that the closed-loop response performance degradation is distributed as a chi-squared random variable when no model bias is present. This result also indicated a direct relationship between the mean and variance in the expected performance degradation for that condition, thereby guaranteeing reliable results if the relationship holds. Experimental results confirmed this relationship in the case of no model bias and suggested a similar result may hold when model bias is present. The D-optimal input signal was realized by a weighted sum of sinusoids. Because the period of one of the sinusoids in the second order case was several times that of the experiment length, this method gave poor results. While we do not generally condemn either alphabetic optimal or input methodologies with sinusoidal time domain realizations, such problems detract from their optimal properties. For the “bias” and “variance” input designs of eqs 20 and 22, realization of the desired input spectrum was achieved by finding an ARMA time series representation with the desired frequency characteristics and then clipping a realization of the series to produce a binary input. This method was found to give adequate results. It is surprising that the issue of obtaining practical time domain realizations of desired input spectra has not received more attention in the system identification literature. Of the designs considered, the bias optimal input (eq 201, denoted M3, gave consistently the best performance for situations involving both correct model structures and model mismatch. Examples of approximate improvements in average controller performance using the M3 design for the cases considered are shown in Table 1. The percentage improvement values shown in this table are the percentage increases in mean(APerf) and Std(APerf) resulting from the use of an input other than M3 (i.e., M1, M2, M4, or M5) over the corresponding values obtained using input M3. It is emphasized that the PRBS switching times were chosen in an ad hoc manner and superior performance
M1 PRBS, T,, = 5 57 (62) 51 (57) 44 (70)
M2 PRBS, T,, = 10 27 (38) 29 (21) 10 (24)
M4 variance (eq 22) 0 (8)
41 (38) 4 (29)
M5 D-optimal 58 (74) 26 (17) 22 (18)
from such inputs could be expected if either this or the switching probability were chosen according to some criterion of optimality as in Tulleken (1990). However, we are not aware of any studies which have used optimal PRBS input signals when model bias is expected. In the face of model mismatch, a properly designed input signal will ensure that the model is matched to the true system a t the relevant frequencies for controller performance. This is clearly the most important factor in ensuring good controller performance. Careless choice of an input test signal can clearly lead to much poorer controller performance and predictability. Better models can usually be identified by increasing the amplitude of the input signal. However, in the case of a biased model with one-step-ahead prediction, higher amplitudes actually led to an increase in the standard deviation of the controller performance degradation criterion when significant frequency mismatch existed. This reemphasizes the need for ensuring proper frequency content in the test signal. The most important limitation of the above results is the dependence of Ljung‘s input designs (eqs 20 and 22) on the disturbance model. Knowledge of the correct noise structure has been assumed in all of our examples and it is not yet clear how sensitive the results are t o deviations from this assumption. Zarrop (1979) has given several examples showing that D-optimal designs are only mildly sensitive to parametric deviations in the disturbance model. Interactions between the choices of input design method, parameter estimation scheme, and controller design also deserve further study.
Acknowledgment Adrian Clark performed many of the Nyquist ellipse computations. The authors gratefully acknowledge the financial support provided by Queen’s University, the Natural Sciences and Engineering Research Council of Canada, and the Shell Development Company.
Nomenclature a = white noise variable A = amplitude of input signal A = coefficient in eq 26 B = coefficient in eq 26 C(o)= frequency weighting function Cov[ 1 = variance-covariance matrix G = set of constraints D = disturbance Det[ ] = determinant e = white noise variable E = statistical expectation f = closed-loop performance degradation signal F(dw)= derivative of closed-loop tracking error g = gain of dynamic system G(z)= process transfer function H(z) = disturbance transfer function J = input design cost function
2666 Ind. Eng. Chem. Res., Vol. 33, No. 11, 1994 k = number of sampling intervals for prediction lead time I( ) = likelihood function lim = limit M = Fisher's Information Matrix n = order of dynamic system N = total number of observations APerf = closed-loop performance degradation measure P = number of parameters 8 T(z) = vector containing G(z) and H ( z ) Tr[ 1 = trace U = manipulated input signal W = weighting matrix Y = measured system output z = unit forward shift operator Greek Letters a = constant chosen to satisfy input power constraint y = transformed transfer function parameter 6 = transfer function parameter 8 = parameter vector , I= eigenvalue p = constant chosen to satisfy input power constraint @ = spectral density function o = frequency Vo = gradient with respect to system parameters Subscripts B = bias contribution c = controller d = disturbance i = imaginary component M = for a specified model N = based on N observations r = real component
= based on an infinite number of observations
Abbreviations PRBS = pseudorandom binary sequence SISO = single input'single output STN = signal-to-noise ratio = with
respect to
Appendix. Optimal Input Spectra and Time Domain Representations A.l. No Model Bias: Bias Optimal Design. The design procedure follows Ljung (1987). The input spectrum is given by eq 20. A controller is to be designed so that the step response of the closed-loop system is the same as that of the open-loop system. We approximate the z transform of a step function by 141 - 0.982-'1. The weighting function based on a pole placement design procedure for this scenario is given by
(A.1) The disturbance spectrum is
Ut = 1.122Ut7,_,- o.470U,-2 - O.653Ut+ 4-a, + 0 . 7 4 1 ~ , - (A.3) ~
A.2. No Model Bias: Variance Optimal Design. Following Ljung, it can be shown that the desired spectrum is the square root of the expression in eq 20. The time domain realization of this spectrum was approximated by a lower order time series model and the parameters determined by a numerical optimization procedure, with the following result:
U, = 1.955Ut-,
-
+
+
1.167Ut7,_, O.1513Ut+ a, - 0.60at-, (A.4)
A.3. No Model Bias: D-Optimal Design. The spectrum of the input requires the numerical optimization of a determinantal equation (Zarrop, 1979). The optimal input consists of a sum of two sinusoids:
+ (0.333)"'
sin(0.0034t)
(A.5)
Literature Cited
Superscripts opt = optimal T = transpose of a matrix - = departure from average value A = estimate - = average = discrepancy between estimated and true values
wrt
The time domain realization of the input spectrum is
Ut = (0.667)"' sin(0.157t)
sp = set point t = at time t V = variance contribution 0 = true value
-*
(A.2)
Astrom, K. J. Introduction to Stochastic Control Theory;Academic Press: New York, NY,1970. Astrom, K. J.; Wittenmark, B. Computer Controlled Systems: Theory and Design; Prentice-Hall Englewood Cliffs, NJ, 1984. Box, G. E. P. Some Theorems on Quadratic Forms Applied in the Study of the Analysis of Variance Problems: Effect of Inequality of Variances in One-way Classification. Ann. Math. Stat. 1954, 25,290-302. Box, G. E. P.; Jenkins, G. W. Time Series Analysis: Forecasting and Control; Holden Day: San Francisco, 1976. Box, G. E. P.; Draper, N. R. Empirical Model Building and Response Surfaces; Wiley: New York, 1987. Davies, W. D. T. System Identification for Self-Adaptive Control; Wiley Interscience: London, 1970. Espie, D.; Macchietto, S. The Optimal Design of Experiments. AZChE J. 1989,2,223-229. Gevers, M. Connecting Identification and Robust Control: A New Challenge. ZFACIZFORS Symp. Zdentif. Syst. Parameter Estim. 1991,9th, 1-10, Goodwin, G. C.; Payne, R. L. Dynamic System Identification: Experiment Design and Data Analysis; Academic Press: New York, NY,1977. Goodwin, G. C.; Gevers, M.; Mayne, D. Q. Bias and Variance Distribution in Transfer Function Estimation. ZFACIZFORS Symp. Zdentif. Syst. Parameter Estim. 1991,9th, 952-957. Goodwin, G. C.; Gevers, M.; Ninness, B. Quantifying the Error in Estimated Transfer Functions with Application to Model Order Selection. IEEE Trans. Autom. Control 1992,37,913-928. Gupta, N. K.; Mehra, R. K.; Hall, W. E. Application of Optimal Input Synthesis to Aircraft Parameter Identification. J. Dyn. Syst., Meas. Control 1976,92,139-145. Kedem, B. Binary Time Series; Marcel Dekker: New York, NY, 1980. Kiefer, J.; Wolfowitz, J. Optimum Designs in Regression Problems. Ann. Math. Stat. 1969,30,271-294. Ljung, L.System Identifiation: Theory for the User;Prentice Hall; Englewood Cliffs, NJ, 1987. MacGregor, J. F.; Harris, T. J . The Exponentially Weighted Moving Variance. J . Qual. Tech. 1993,25,106-118. Mareels, I. M. Y.; Bitmead, R. R.; Gevers, M.; Johnson, C. R.; Kosut, R. L.; Poubelle, M. A. How Exciting Can a Signal Really Be? Syst. Control Lett. 1987,8, 197-204. Mehra, R. K. Optimal Inputs for Parameter Estimation in Dynamic Systems-A Survey and New Results. IEEE Trans. Autom. Control 1974,AC-19, 753-768.
Ind. Eng. Chem. Res., Vol. 33,No.11, 1994 2667 Rippin, D. W. T. Statistical Methods for Experimental Planning in Chemical Engineering. Comput. Chem. Eng. 1988,12,109116. Rivera, D. E.Monitoring Tools for PRBS Testing in Closed-Loop System Identification. AIChE National Meeting, Miami-Beach, 1992,Paper 131-d. Rivera, D. E.; et al. An Industrial Perspective on Control-Relevant Identification. Proc. ACC 1990,Paper FA11, 2406-2410. Rousseeuw, P. J.;Leroy, A. M. Robust Regression and Outlier Detection; Wiley: New York, NY,1987. Soderstrom, T.; Stoica, P. G. Instrumental Variable Methods for System Identification; Springer-Verlag: Berlin, 1983. Tulleken, H. J.A. F. Generalized Binary Noise Test Signal Concept for Improved Identification-Experiment Design. Automatica 1990,26,37-49. Tulleken, H. J. A. F. Grey Box Modeling and Identification Using Physical Knowledge and Bayesian Techniques. Automatica 1993,29,285-308. Wahlberg, B.; Ljung, L. Hard Frequency Domain Model Error Bounds from Least Squares-Like Identification Techniques. IEEE Trans. Autom. Control 1992,37,900-912.
Watts, D. G.; Minich, G. M. Design of Experiments for Estimating Process Dynamics. Preprint No.1971-52,Department of Mathematics and Statistics, Queen's University, Kingston, Ontario, Canada, 1971. Wilkinson, R.; Harris, T. J.; Davis, J. The Numerical Inversion of the Characteristic Equation with Applications to Positive Quadratic Forms in Normal Variables. Submitted for publication in Commun. Stat., 1993. Wittle, P.Hypothesis Testing in Time Series Analysis; Thesis, Uppsala University, Almqvist and Wiksell: Uppsala; Hafner: New York, NY,1951. Zarrop, M. B. Optimal Experiment Design for Dynamic System Identification; Springer-Verlag: Berlin, 1979. Received for review January 11, 1994 Accepted July 7, 1994 @
Abstract published in Advance ACS Abstracts, October 1, 1994. @
