- Email: [email protected]
Integrated identification and robust control
J. Proc. Cont. Vol. N, Nos. 5 6, pp. 431 44(I, 1998
ELSEVIER
1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 095t)-1524;98 $ see front matter
PII: S0959-1524(98)00027- 4
Integrated identification and robust control Brian L. C o o l e y and Jay H. Lee* Department of Chemical Engineering, Auburn University. Auburn. AL 36849-512 7, USA
A framework for integrated identification and control is presented. As part of this framework, frequency domain uncertainty bounds are derived for robust stability tests, a robust stability test for elliptical bounds is developed for SISO systems, a methodology for estimating controller performance is derived, and an optimal experiment design methodology for control-relevant identification is outlined. An example is presented to illustrate how the tools of the framework fit together. ~(' 1998 Elsevier Science Ltd. All rights reserved Keywords: identification for control; model-based control; robust stability; optimal experiment design
In this paper, the need for an integrated identification and robust control methodology is addressed. In recent work, this phrase has been used to refer to supplying frequency domain uncertainty bounds for use in robust controller design and analysis. Here, the phrase is used to refer to a model-based controller design methodology in which the various stages of identification are conducted with respect to the final goal of high performance control. The following aspects form the core of the methodology, which is depicted in Figure 1.
arriving at an uncertainty description and using uncertainty information for designing test signals. The first stage of the proposed methodology is to translate the system's statistical information and any available data into a model and an associated uncertainty description. Here, what we arrive at are ellipsoidal bounds on the frequency response of the identified model. While this uncertainty description may be approximated for use in existing closed-loop analysis and robust controller synthesis tools such as/~-theory, the use of such approximations may at times lead to overly conservative results, as we will demonstrate with an example. In such cases, it is necessary to develop an alternative method for testing robustness and designing robust controllers. A straightforward method for directly using the developed ellipsoidal bounds to test robust stability of SISO systems is presented. The complexity of the robust stability criterion for MIMO systems makes extension of the direct test of the elliptical bounds to M I M O systems difficult, and such extension is left for future work. The second stage in the proposed methodology is to use the identified model and uncertainty description to design future experiments. The focus here is to design control-relevant experiments, i.e. experiments that generate data most useful for developing a high-performance model-based controller. We adopt a statistical framework in which future experiments are designed based on past information. The test signal is set equal to the argument that minimizes the conditional expectation of a closed-loop performance index. This paper is arranged as follows. In the following section, the problem of deriving frequency domain bounds is addressed. 'Soft' bounds are derived in a manner similar to that of previous studies 12. In 'Robust stability tests', the problem of testing robust stability
•
Estimating uncertainty. Since model error causes controller performance to deteriorate, it is important to have a useful uncertainty description for predicting stability and performance. For example, frequency domain bounds on the system frequency response would be useful in /z-synthesis and /zanalysis. • Shaping uncertainty. It is intuitive that the test signal used in identification affects the model error distribution and, hence, the controller performance. Ideally, the test signal should be used to influence the error distribution in an optimal (for control) way. • Accounting for uncertainty. Once an uncertainty description is provided, it should be used to analyse controller performance and possibly to improve controller performance. The problem of integrated identification and control covers a broad range of disciplines, and by no means does this paper represent a complete solution. In this paper we will address the two most important stages in an integrated identification and control methodology: *To w h o m c o r r e s p o n d e n c e s h o u l d be a d d r e s s e d . Fax: (205)844-2063; e-mail: jhl!a e n g . a u b u r n . e d u
431
Integrated identification and robust control." B. L. Cooley, J. 14. Lee
432
__~
,,st
Uncertainty Quantifier
Identification Algorithm
YES P
Analysis
NO Control-Relevant Test Input Generator Figure 1 Integrated identification and control methodology
using the derived bounds is approached. It is shown that the derived bounds can be directly tested for singleinput, single-output (SIS®) systems, but only conservative tests are available for multiple-input, multipleoutput (MIMO) systems. In 'Performance tests', the problem of estimating performance based on parameter estimates is addressed. In 'Experiment design', a method for using prior information to design the test signal in a control-relevant manner is developed. Finally, 'Numerical example' presents an example illustrating how the pieces of the proposed integrated identification and control methodology fit together.
Estimating uncertainty Deriving time-domain bounds The underlying system is assumed to be described by a finite expansion in a set of orthonormal basis functions
fk(q):
fk(z) = ~ z
(1 - az) k ~ f a ~ i [a[ < 1
(5)
Other basis sets exist, but the two mentioned above are perhaps the most widely used. Note that model (1) makes sense only for stable systems. Remark: The assumptions on v(t) are not as restrictive as they may seem. If v(t) can be described as v(t) = H(q)e(t), where e(t) is a normally distributed i.i.d, random sequence of zero-mean and covariance I,y with H(q) known or estimated, the model can be put into the form of (1) by filtering the output with H-~(q). This is known as 'prewhitening'. The estimated model then has the form H-l(q)G(q), where G(q) is the system transfer function. Given (1), an estimate of ® can be obtained via least squares from a set of N input-output data. N
(9 = P Z
¢(t)yr(t)
(6)
t=l n
y(t) = ZOgfk(q)u(t) + v(t) = ore(t) + v(t)
(1)
-1
k=l
(7)
®r = [01 ... 0~]
(2)
q~(t) = [fi (q)ur(t)...f~ (q)ur(t)] 7"
(3)
In the above, u E R n" is the system input vector, y E R'y is the system output vector, Ok E R ny×n" are the (realvalued) expansion coefficients, and v C R "~ represents the effects of noise and disturbances. The sequence {v(t)} is assumed to be a Gaussian independently and identically distributed (i.i.d.) random sequence of zeromean and identity covariance. The bases fk(q) in the forward shift operator q (i.e. q-lu(t) = u(t - 1)) have the following orthonormality property:
fk(d'°)fi(e-J°~)dw =
0 k ¢ l
(4)
--Tr
Examples of such bases include fk(z)=z -l impulse response models) and Laguerre bases:
(finite
If v(t) ts uncorrelated with ~(t), then O will be an unbiased estimate of ® (i.e. E® = ®). Furthermore, t~ has the following covariance: ng blocks
Is) O" In the above, vec(.) is an operator that stacks the columns of a matrix in a vector:
vecIa'a211=[a'
Integrated identification and robust control.'B. L. Cooley, J. H. Lee Deriving time-domain bounds
Then for w = Ax C R"", nw<_nx, A c R ""×'' with full rank,
'Soft' uncertainty bounds may be derived using the X2 distribution: bt = vec(O - o)rE(elvec((2) -- O) --~ X2Co)
(10)
where p = n ny n, is the number of parameters. Equation (10) simply states that bt has X2 distribution with p degrees of freedom. This fact enables the construction of confidence bounds:
P[bt < X2×(p)] = y
433
wT(AZAT)-Iw
T T T gkl(ATSkiZoSklAT)-lgkl < X2y(p)
= FRe(dk'-Gk') ] gkl
Deriving frequeno,-domain bounds
(16)
1
Proof" See Wahlberg and Ljung 4. Using Lemma 1 with ( l l ) and (14), frequency response bounds can be derived as follows:
(11)
In the above, 7 is the probability (or confidence) that the true parameters lie within the specified ellipsoidal set, and X2(p), which depends on 7, determines the size of the set. For example, if 7 = 0.95, then X~(2)= 5.99 (see Papoulis3). Equation (11) provides a simple method for bounding the time domain parameters.
<
LIm(Gkt
(17)
(18)
Gkt) J
The above inequality describes confidence ellipses in the Nyquist plane. The ellipses are centered at Gkt and have axes defined by A TSkl Z (0Sklr A rT"
Remarks." 1. The bound (17) is based on the joint distribution of the model parameters. In other work 12, bounds are drawn based on the marginal distribution of
gr:
To test the stability of a particular controller, it is convenient to translate the time domain bounds into the frequency domain. Notice the following
Gkt(e/', (0) = A(d°')S~tvec(O)
(121)
where Gkl(d% O) is the k,l element of the frequency response matrix, and Skz is a selector matrix that selects the appropriate elements of 0 to construct the k,l element. A,(e #°) has the following form:
A(d '°) ~ [A, (d'°)... A.(d°')]
(13)
where Ak are the frequency responses of the corresponding orthonormal basis function. For example, if fk(q) = q k then Ak(# '°) = e -jk'°. Notation: In the following the arguments #~o and 0 will be dropped when it will not cause confusion. The frequency response of the underlying system will be referred to as G = G~d', ®), and that of the estimated model as d = G(d °', 0). It is useful to split the frequency response into real and imaginary parts:
r-e/<,l = [ Im(Gkt) J
Im(A)
Sklvec(O)
(14)
AT
1" A T T gkl( TSklEOSktAr)
-1
~ gkt <_ X•(2)
(19)
The result (17) suggests that the above bound actually has confidence V' such that X~/(P) = X2(2). If 7 is significantly less than 1 (e.g. V= 0.95) or n is large, the set defined by (19) can be considerably smaller than required by (17). 2. Using the bounds in (17) causes the loss of information about the correlation with respect to both frequency and the different elements of G. In the context of a robust stability test, losing the frequency correlation introduces no conservatism since the test need fail at only one frequency. However, losing the correlation between elements can lead to very conservative results since the multivariable Nyquist criterion will strongly depend on the structure of the uncertainty between different elements of the frequency response matrix. While robust stability of SISO systems can be tested using the bounds in (17), robust stability of MIMO systems depends on all the elements at a particular frequency. Using (12) and Lemma 1, it is straightforward to show:
gT(sN(.)ST)-'g < X2×(p)
(20)
where where Re(.) extracts the real part of a complex number and Im(.) extracts the imaginary part. To translate the time domain bounds to frequency domain bounds, the following result 4 is needed: Lemma 1: Let x E R "~, Z > 0 C R n*xn." and xrE Ix< 1
(151)
I Re(Gll - G l l ) Im(Gll. - Gll) g
Re(A)&l Im(A)Sli
(21)
S = [ I m ( G .... -- a < . . )
Re(A)Sn,n. Im(A)Sn,..
434
Integrated identification and robust control. B. L. Cooley, J. H. Lee
Element-by-element circular bounding.
The bounds in (20) (or even (17)) are difficult to use for M I M O systems. For use with standard analysis and synthesis tools, the bounds must be approximated with circles. It is instructive to examine the singular value decomposition (SVD) of A Skiv A r. r"
As an alternative, it is possible to retain information about the structure of the uncertainty in the frequency response elements for a given output by noticing that: AGrowk = ( ( G -
TSkl]~0
TTT : !U1.~ f ffl,k'o f2,k'O] u~TI ATSkI~aoS~IA Ukl
~
(22)
J
(23)
Idk, - Gktl <_ Cx~Co)m,~, Define
I Ac'kl 1 " A
1 E~/2(AGrowk)T
(24)
L &.,k..J
In the above, (A)rowk is the kth row of A (see (25)), and ® represents the Kronecker product• Note that A (not At) appears in (27). The following bound may be derived:
(28) In the above, AG is the complex conjugate of AG. This bound, although it incorporates information about the uncertainty structure, is also conservative. This is easy to see since, for a SISO system, the bound approximates an ellipse (17) with a disk ((28) is a magnitude bound on a complex scalar for SISO systems). Unfortunately, the disk is larger than the bound defined by (23). However, (28) contains information about the relationship between the frequency response elements for a given output. Consider the following singular value decomposition: (a ®
AGrowk ~ [Gkl -- Gk,... dk,~. - Gkn,]
O'l,kl
0
-"
.
"
(27)
AGrowk((A ®ln.)P(A* ® I,~.)) l(AGrowk)T < X2v(p)
Gkl
where al,kl ~ Cr2,kl. Thus, lies in an ellipse of center (~kl with axes ul and u2, of lengths ~ and ~v/~,~/. The smallest circle that overbounds this set has radius (see Figure 2). This circle is described by
Ac,k =
G)T)rowk(AT @In.)
(25)
0
I,,.)P(A* ® i,,.)
(26)
(29)
where * denotes the complex conjugate transpose, Us is unitary, and E~ is diagonal. Define
As,k =
~ac, k
U, sG
=
I As'kl " ]
A
L A,'k.°
l
~l/2G(AGrowk) T
(30)
"
0
0
-..
0
(71,kn.
It is easy to show via (23) that
]Ac,kzl
< 1, l = 1. . . . . n,.
Row-by-row structured bounds. Obviously, approximating the bounds of (17) with (23) is conservative when ~rl,kl > ~rZ,kt since the circular bounds include values for Gk/ not in the elliptical bound• However, a greater criticism of using these bounds is that they are independent bounds on each element of the frequency response matrix and therefore do not incorporate information about the structure of the uncertainty. This is potentially very conservative. Im
In the above, As,kt, I = 1. . . . . n, are complex random variables. It is easy to show via Lemma 1 that (28) implies that IAs,kll _< 1, l = 1. . . . . n, Remarks. • Both (23) and (28) are conservative bounds. However, the relative conservativeness of these bounds depends on the characteristics of Eo. If this matrix is such that there is little correlation between elements of the frequency response matrix, then it is expected that (23) will be less conservative than (28). On the other hand, when such correlation exists, (23) can be very conservative in comparison to (28). In general, both bounds should be computed, and the structured singular value should be taken to be the minimum value for the two bounds.
O-lU 1
Robust stability tests SISO systems
Re Figure 2 Overbound approximation to (17)
While it would be simple to use /z-theory (i.e. existing tools) by approximating the elliptical bounds in (17) with over-bounding circles, such an approach would clearly lead to conservative results. Instead, we will
Integrated identification and robust control." B. L. Cooley, J. 14. Lee derive a robust stability test that explicitly uses the elliptical bounds by translating the stability condition for SISO systems into a test based on the elliptical bounds. For SISO systems, the robust stability criterion simply asks whether G C = - I for any model in the model set and a given controller C. This simple criterion is easy to test upon considering the following representation of the product of GC: G C -- Re(C)Re(G) - Im(C)Im(G)
(31) + (Im(C)Re(G) + R e ( C ) I m ( G ) ) j
!Re(c) Im(C)
1 rRe(oc)l LIm(GC)]
Im(G)
Re(C)
=
(32)
Note that the (frequency-dependent) transformation matrix CT maps from the complex plane to the complex plane. Define
g'
A [Re((G-G)C)] O)C)
(33)
435
ellipsoidal set centered at G C with axes determined by (CTATE(.)ATCTT) that does not contain G C = - I is determined by T
gr
T T (CTAT]~(')ATCT)
I
gc
<
rain h
(39)
Remarks. •
1. Determining the maximum confidence is very desirable since Y is usually chosen somewhat arbitrarily. Knowledge of the maximum Y would allow the engineer to make an intelligent choice between accepting the model or collecting more data to refine the model. For example, if the acceptable confidence is 7=0.95, then the engineer might decide not to retest if it were determined that min~ b = X~.94(F/). Ideally, the y associated with minoj b could be monitored on-line and used to decide when to stop testing. 2. The robust stability test outlined above is similar in concept to that of a number of previous papers 5 7. Note that the above robust stability test and the tests described in the cited works are exact tests (i.e. no approximation of the bounds is required). M I M O systems
The mapped bounds (n,, = ny = l) become gc Zx, g' <- x?(n)
A
~g, =
(34)
T
T
(35)
CTA T~(.)A TCT
To check the bounds for robust stability, define = L
Im((~C)
z_~ g'
Re(dC) + l Im((~C)
(36)
Since and (~ and C are known, b is a known function. Robust stability can be determined by checking the following inequality at each frequency: 9
b < X~,(n)
(37)
If (37) holds at any frequency, then the uncertainty set defined by (34) contains the point G C = 1 in its Nyquist band, and robust stability does not hold with confidence % The above stability test can be extended to determine the maximum confidence for robust stability. Notice that b >_ m i n h
(38)
2 Thus, if y is selected such that X× < min~ b, then (37) will fail at every frequency. In addition, the largest
To avoid conservativeness, MIMO systems should be tested in a manner similar to that outlined in the previous section for SISO systems. Unfortunately, the multivariable Nyquist criterion, which asks whether the map of det(& + GC) contains 0 for any model within the model set characterized by the bounds in (20), is complex, and it is difficult to extend the simple technique developed for SISO systems to MIMO systems. As an alternative, structured singular value (la) theory can be used to test stability./~-Theory is a powerful tool for checking robust stability, but the test assumes that the frequency domain bounds are disk-shaped. While the bounds of (17) cannot be used directly, those of (23) or (28) can be used. Since the relative conservativeness of the two bounds depends on the structure of Z(.}, a stability test should be performed for each type o[" bound, and the least conservative result (smallest values of the structured singular value) should be used. The bounds in (23) can be used as lbllows. First, uncertainty descriptions can be constructed for each output. It is apparent from (24) that
AGro,~k -- [1-
0
.11][{£ ." ,6,j,i
•.
0 ]
•
0
0
A,. k.,,
(4o)
Integrated identification and robust controL B. L. Cooley, J. H. Lee
436
The uncertainty descriptions for each output can be combined as shown in Figure 3 with Z = Zc (see below) and D = I,,,y. This system can easily be put in standard M - A form 8 with
M = Wl,sC(Iny -~- G C ) - I W2
,,,',,
M = Wl,cC(In~ +
2
(~C)-I W2
"½
f/-}
=
,,, blocks
with the additional definitions
ny blocks; ",
Us zx
".°
"..
"•
"..
(42)
il
"..
• "•
1!;
•,
"•
•
"..
". •
0
...
0
A
0z~!0
n. elements
0
''"
,.
*,
••"
0
...
0!
(43)
0
(44) Ec.,~
For the bounds of (28), the uncertainty description for a particular output Yk can be written as follows: Akl
0
...
0 ]
..
".
".
".
0
0
Akm
Aarowk = [ 1 . . - 1 ] 0
x
•
" "
1
(45)
V/ 2(p)Zs1/2 UsT
The uncertainty each output and = Zs and U = in standard M - A
descriptions can be constructed for combined as shown in Figure 3 with Us (see below). This can easily be put form with
.
The uncertainty description for the M I M O system
0 ]
~
1
(49)
A structured singular value test can be conducted for M of the form in (41) and for M of the form (46). Since both tests are based on conservative bounds, the smallest structured singular value obtained from the two tests can be used as an upperbound on the structured singular value that would result from testing the true ellipsoidal bounds• R e m a r k : In 'Numerical example' an example will be presented to demonstrate the potential conservatism of the robust stability tests based on (23)/(41) and (27)/ (46). It is important to note that such conservatism arises from force-fitting the ellipsoidal bounds into bounds that can be used in standard/z-analysis. There are two ways to reduce (or even eliminate) the conservatism• First, less conservative approximations of the ellipsoidal bounds may be sought. Second, alternative robustness tests that directly use the ellipsoidal bounds may be developed. It is our belief that the latter approach is the correct one, and it is the approach we used for SISO systems (see 'SISO systems')• However, for SISO systems, the stability test involves only checking whether a point lies within an ellipse• For MIMO systems, the stability criterion is much more complicated (det(I,~ + G C ) = 0) and not so easily checked. Clearly, the search for such alternative robustness tests deserves a paper in its own right and is beyond the scope of the current paper.
Performance
Figure 3
...
n~ blocks
n. elements
0
(48)
0
0
$ ----
y
]~e,l
.•
Z~
with the additional definitions ...
i1 .•
ny b]Cocks
ny blocks
iS[1
(47)
(41)
I"]} •
(46)
tests
While robust stability is desirable, the ultimate goal is good performance. For example one might be interested in obtaining good mean-squared error performance: J = EeT(t)e(t) = tr(E{e(t)er(t)})
(50)
Integrated identification and robust controk B. L. Cooley, J. H. Lee In the above, e represents the difference between the reference signal and the system output (see Figure 4). It is straightforward to show
437
An interesting result is that minimization of the righthand side of (54) is equivalent to minimizing the weighted trace of the parameter covariance matrix: min J ~ min W'Zo
e=S(Iny + (G - (~)(}-I/})-1 (r - d)
(56)
(51) A formula for the weight W can be found elsewhere 9.
Remarks." In the above S = (I.y + GC) -I and H = I . , - S. d is assumed invertible for convenience. Using the above expression and Parseval's relationship, it is easy to show that (50) becomes J = t r ( I~
f~ E{S(I< +(G-O)G
1/})-1 (52)
1. The controller does not explicitly aflpear in (54). Instead, the sensitivity function S is specified. Since the model G is available, determining S for a given C is straightforward= Note that it is a common practice to specify S and translate it to C once G is obtained. 2. Care should be taken in using (54) since it only holds when the closed-loop system is stable.
x q%(/< + (G - G)d-l/-})-*S*}dw)
Experiment design In the above, the superscript * represents complex conjugate transpose (A-* is the inverse of A) and qbn is the spectrum of rl. Equation (52) directly measures the expected performance of the controller. However, it is difficult to evaluate because the unknown plant G enters into the expectation in a nonlinear fashion. When the plant estimate is close to the true plant, (51) may be linearized about (~, which is the best estimate of the plant: e
+
(53)
-
Thus, for unbiased (~
(54)
The stability test and the performance test described in the preceding sections provide criteria for accepting or rejecting a model. When the model is rejected, more data must be collected to refine the model. The information gained from previous experiments can be used to extract more control-relevant information from the system through test signal design. It is intuitive that the optimal experiment design will depend upon the system dynamics. Thus, it is necessary to employ an iterative scheme where the design is computed based on a prior distribution of the plant parameters. The test signal design problem is formulated as follows. First, assume that a prior distribution on the plant is given in the form of a nominal model (2)k, and its covariance E(.)k. This prior might be specified by the engineer, or it might come from the estimates provided from a previous experiment. When new data are collected, the model can be updated as follows:
N
) 1 (57)
The second term in the integrand of (54) is a linear combination of E{((~ - G)*(G - G)} which is given by Ok+l = Z (Ok -1 + I E{(G - G)*(G - G)} =
(P(t)yr(t)+ E6~It2)k
(581
ny(A r ® I..)*P(A r ® I..)
(55) Notation: Here, the notation is extended to incorporate In the above ® represents the Kronecker product.
r
_
~
the iterative design by defining, for example, G~-G(d°~,~)k). This is not to be confused with GrowkThe test signal is designed to minimize the conditional expectation of the performance index in (52) subject to any constraint on the inputs:
y u~kP+tl= arg min
Figure 4 The closed-loop system resulting from identification of model G and constructionof model-basedcontrollerC
Uk+l E U
E{JlOk, I2~.~}
(59)
(60)
Integrated identification and robust controk B. L. Cooley, J. H. Lee
438
In the above, Ug+l is the test signal used in the (k+ 1)st experiment. U is a set that expresses constraints on the test signal and is usually characterized by magnitude and rate bounds. Note that J depends on the plant, the reference signal/disturbance spectrum, and the model obtained after the experiment: J
=
tr
Sk+l(/., + (G
x O,(I,> + (G
-
X S~+l](~kZOk
-
--
^
^-I ^
Gk+l)Gk+lHk+l
)-I
Gk+l)Gkll/~k+l) -*
})
doo
(61) The plant and the model are unknown and treated as random variables. The minimization (59) is difficult to solve because the objective function contains the unknown plant and the model, which enter in a nonlinear fashion. It is convenient to linearize the error expression, this time with respect to both the plant and the model. The linearization is taken around the best estimate of the plant, Gk: e ~ (Sk+l + Sk+l((}k+l - G)G/1/Qk+I)~
(62)
Note that, in general, Sk+l and /}k+l depend on (~k+l. However, if a direct synthesis approach to designing the controller is taken (i.e. the sensitivity function is specified), then Sk+l and Hk+l become design variables (i.e. their values are specified by the engineer). The minimization in (59) becomes U~kpt' = arg min tr
~
E Sk+l (dk+l - G) -Tr
^
1
^
^
*
^
x G k g k + l ~ ( G k I-I~k+l
(63)
-- Gk+l) S~+I]OkZO k d o
It is straightforward to show that the above objective function is simply a linear combination of the elements of E(() - G)*(G - G), and can thus be related to the test signal through (55). As a result, the above optimization can be solved, but the potentially large dimension of Uk makes computation an issue. A suboptimal way to compute this optimization is to divide the test signal into small portions and solve for each portion. Furthermore, (63) is equivalent to minimizing the weighted trace of the parameter covariance matrix: min J ~ min WZo
resulting model (and the uncertainty) by denoting it with an index corresponding to the experiment being performed (i.e. &+l). Inclusion of uncertainty in the controller design greatly complicates the experiment design, however. Note that, in this case, Sk+l is not known prior to the experiment and hence is a random quantity. To solve the problem rigorously, one must embed into the expectation the dependence of S upon the future model and the uncertainty. This topic is worthy of a paper in its own right, and we defer an indepth study to future research. In practice, since the objective is only^to improve an experiment design, assuming Sk+l = Sk will generally be sufficient.
(64)
A formula for the weight W can be found elsewhere 9. Remark." One could potentially consider a more complicated robust control system design method, rather than a direct synthesis method (for which the closedloop functions are fixed). We implied that the desirable closed-loop nominal response may depend upon the
Numerical example Preliminaries In this section, we will test the integrated identification and control methodology discussed in this paper on a nonlinear high purity binary distillation column model. The model for the column has been set forth previously l° and has been considered in a number of recent works 11,j2. Briefly described, the column consists of 25 trays, a reboiler, and a total condenser and is designed to separate a binary mixture of methanol and ethanol into product streams of 99% purity. The two controlled variables are the temperatures on trays 21 (T21) and 7 (TT). The two manipulated variables are the reflux flow (L) and vapor boilup (V). The nominal values are T21=614.1°R, T7=636.8°R, L-3.3835mol/min, and V = 3.856 mol/min. The column model includes tray-bytray energy balances and tray hydraulics. The objective in this example is to use the integrated identification and control methodology to identify models for use in building controllers. The controllers are specified by designing a desired closed-loop response. Here, the complementary sensitivity function is fixed so that the nominal response of the outputs to the reference signal is first-order with time constant 4: /} = (1 - e-1/4)z -1 1 e_l/4z_ 1 12
(65)
This corresponds to the following controller: C(z) - (1 - e-1/4)z -' ~ 1 1 -
z -1
(z)
(66)
The above controller is used when the model contains no zeros outside the unit disk. When the model contains zeros outside the unit disk, the model is separated into all-pass and a minimum phase portions as follows: G(z) = GA (Z)GM(Z)
(67)
where Ga(z) is all-pass, and GM(Z) is minimum phase. See Morari and Zafiriou 8 for details on this inner-outer factorization. The following controller is used:
Integrated identification and robust control." B. L. Cooley, J. H. Lee
C(z) = (1
-
15
e 1/4)2-1
i --- ~
439
GMi (z)
'l
(681) 10
While the above controller may give a different nominal sensitivity function from the one used in the experiment design, the sensitivity functions are usually simillar. The above controller will maintain closed-loop stability when the model contains unstable zeros, however. Note that (66) and (68) are controllers with integral action. To obtain the models, preliminary experiments are performed to determine the settling time (around 240 min). Based on these results, we decided to sample the column every 4min and model it using a finite impulse response model consisting of 60 FIR coefficient matrices. Two sets of experiments are performed on the column. In the first set, 1000 samples are generated using an independent random binary signal (RBS) for each input channel of deviation +0.012mol/min from the nominal values. In the second set of experiments, the experiment design methodology of 'Experimental design' is used. The experiments are designed by breaking each experiment into four sections, each consisting of 250 samples. The first 250 samples are generated using an independent RBS for each input of deviation +0.012mol/min from the nominal values. A preliminary model is then identified based on this data, and this model is used to design the input signal for the next 250 samples with each input constrained to a maximum deviation of +0.012mol/min. The model is then updated recursively, and the experiment design procedure is repeated using the new model. A total of 1000 samples are collected (250 RBS, 750 designed). For each experiment, the outputs are corrupted with a zero-mean independently and identically distributed random sequence of covariance (0.1°R)212. Controllers are designed based on the resulting models as described above.
Analysis of stabilio, Uncertainty bounds can be constructed for the identified models as described in "Estimating uncertainty' and used to analyse the stability of the closed-loop system for a given controller. Plots of the structured singular value (/x) for the two uncertainty tests are given in Figure 5. It is interesting to note that (23) provides a tighter /~ than (28) for the RBS trials, while (28) provides a tighter/x than (23) for the optimal design. This suggests that it is crucial to capture the structure of the uncertainty when computing the stability measure. One particularly notable aspect of this trial is that the /x-analysis result seems to indicate that the closed-loop system should frequently be unstable, both for the RBS experiments and for the optimal designs. Judging from Tabh, 1, this may be true for the RBS designs, but it is not true for the optimal designs. Moreover, it is difficult to discern from the figure which design is better, despite the fact that the simulation results proved the optimal design to result in a stabilizing controller more often.
iI
1
t~\//"A'Xv/'41, 0
i'L
fi ~ ,~, , ,A'
~' ~
~,
I l, •
0!f
012
013
014 frequency
0'5
016
0L7
, 01
012
0~3
0i4 frequency
015
i 06
i 07
25
>15
(b)
10 5
O0 25 Ii
2o 1: > 15
!o
L,:
/t
,'
" "
1
5 ji.L/~tl_t'~ll ,'/i~izt/~ x lit k' 11 ,1' , 01
02
03
,i /
04 frequency
,.
05
!'
^ 06
11~ , ,I 07
Figure 5 typical result: (a) structured singular value (ix) for RBS design (solid line, computed from (23)/(4]): dashed line, computed from (2?)/(46)); (b) ix for optimal design (solid line, computed from (41); dashed ]me, computed from (46)); (c) comparison of the RBS and optimal design (solid line, smallest IX for RBS design: dashed line, smallest ls for optimal design)
This indicates the conservativeness of the robust stability tests described in 'MIMO systems'. Clearly, these results motivate the search for an alternative robust stability test that can use the derived ellipsoidal bounds without approximation. One approach for reducing the conservatism in MIMO robustness tests is that of highly structured stability margins, developed previously 13"14 While the concept is applicable to the problem developed in this work, algorithms that explicitly deal with ellipsoidal bounds are not developed in these studies. The development of such algorithms is beyond the scope of this paper and is left lk~r future work.
Comparison q/ RBS and optimal design It is of interest to compare the performance of controllers designed based on the models identified from each experiment. Table 1 shows the results of the two sets of experiments. Clearly, the experiment design technique greatly improved the closed-loop behavior, both by yielding more stabilizing controllers and better performing controllers. Note that the unstable runs are such that the matrix G(0)(~ 1(0) has eigenvalues with negative real part (G(0)) is the steady state gain matrix Table I
Simulation results
No. unstable responses Average ISE (stable)
RBS design
Oplimal design
45 {).7827
8 0.3858
Integrated identification and robust control." B. L. Cooley, J. H. Lee
440 614~
and performance. If it is determined that the model must be refined, the information gained from previous experiments can be used to help design future experiments so that more control-relevant information is collected.
6143E 614
~ (a) ~
614 25 814.2
~'614 15 614 1
References
614,05 614 0
I 100
J 200
I 300
I 400 time
I 500
L 300
L 700
J 800
500
600
700
800
837¸1
307 '"
"~"~~
-
~
~ 3 0 6 g5 {b) r,. 636g >, 636 85 636 8 638 75
0
1 O0
200
300
400 IJrne
Figure 6 (a) T21 output response (solid line, RBS design; dashed line, optimal design; dotted line, setpoint); (b) T7 output response (solid line, RBS design; dashed line, optimal design; dotted line, setpoint)
based on the linearized model). In such cases there is no controller (based on the model) with integral action that can stabilize the closed-loop system 8. The results from a particular trial are shown in Figure 6.
Conclusions This paper develops an integrated identification and control scheme that provides not only frequency domain uncertainty bounds but also tests for stability
1. Goberdhansingh, E., Wang, L. and Cluett, W. R., Robust frequency domain identification. Chem. Eng. Sci., 1992, 47(8), 19891999. 2. Goodwin, G. C., Gevers, M. and Ninness, B., Quantifying the error in estimated transfer functions with application to model order selection. IEEE Trans. Automatic Control, 1992, 37(7), 913928. 3. Papoulis, A., Probabilty, Random Variables, and Stochastic Processes. McGraw-Hill, New York, 1991. 4. Wahlberg, B. and Ljung, L., Hard frequency domain model error bounds from least-squares like identification techniques. IEEE Trans. Automatic Control, 1992, 37(7), 90~912. 5. Tsypkin, Y.-Z. and Polyak, B. T., Frequency domain criteria for LP-robust stability of continuous linear systems. IEEE Trans. Automatic Control, 1991, 36, 1464-1469. 6. Blernacki, R. M., Hwang, H. and Bliattacharya, S. P., Robust stability with structured real parameter perturbations. IEEE Trans. Automatic Control, 1987, 32(6), 495 506. 7. Guzzella, L., Crisalle, O. D., Kraus, F. J. and Bonvin, D., Necessary and sufficient conditions for the robust stabilizing control of linear plants with ellipsoidal parametric uncertainties. In Proc. 30th CDC, Brighton, UK, 1991, pp. 2948-2953. 8. Morari, M. and Zafiriou, E., Robust Process Control. PrenticeHall, Englewood Cliffs, NJ, 1989. 9. Cooley, B. L. and Lee, J. H., Control-relevant experiment design for multivariable systems. Automatica, submitted. 10. Weischedel, K. and McAvoy, T., Feasibility of decoupling in conventionally controlled distillation columns. Ind. Eng. Chem. Fundamen., 1980, 19, 379-384. 11. Li, W. and Lee, J. H., Control relevant identification of ill-conditioned processes. Comput. Chem. Eng., 1996, 20, 1023-1042. 12. Chien, I.-L. and Ogunnalke, B., Modeling and control of highpurity distillation columns. In AIChE Annual Meeting, paper 2a. Miami, 1992. 13. Figueroa, J. L., Desages, A. C., Romagnoli, J. A. and Palazoglu, A., Highly structured stability margins for process control systems. Comput. Chem. Eng., 1991, 15(7), 493-502. 14. Rotstein, H., Desages, A. and Romagnoli, J., Calculation of highly structured stability margins. Int. J. Control, 1989, 49(3), 1079-1092.
ELSEVIER
1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 095t)-1524;98 $ see front matter
PII: S0959-1524(98)00027- 4
Integrated identification and robust control Brian L. C o o l e y and Jay H. Lee* Department of Chemical Engineering, Auburn University. Auburn. AL 36849-512 7, USA
A framework for integrated identification and control is presented. As part of this framework, frequency domain uncertainty bounds are derived for robust stability tests, a robust stability test for elliptical bounds is developed for SISO systems, a methodology for estimating controller performance is derived, and an optimal experiment design methodology for control-relevant identification is outlined. An example is presented to illustrate how the tools of the framework fit together. ~(' 1998 Elsevier Science Ltd. All rights reserved Keywords: identification for control; model-based control; robust stability; optimal experiment design
In this paper, the need for an integrated identification and robust control methodology is addressed. In recent work, this phrase has been used to refer to supplying frequency domain uncertainty bounds for use in robust controller design and analysis. Here, the phrase is used to refer to a model-based controller design methodology in which the various stages of identification are conducted with respect to the final goal of high performance control. The following aspects form the core of the methodology, which is depicted in Figure 1.
arriving at an uncertainty description and using uncertainty information for designing test signals. The first stage of the proposed methodology is to translate the system's statistical information and any available data into a model and an associated uncertainty description. Here, what we arrive at are ellipsoidal bounds on the frequency response of the identified model. While this uncertainty description may be approximated for use in existing closed-loop analysis and robust controller synthesis tools such as/~-theory, the use of such approximations may at times lead to overly conservative results, as we will demonstrate with an example. In such cases, it is necessary to develop an alternative method for testing robustness and designing robust controllers. A straightforward method for directly using the developed ellipsoidal bounds to test robust stability of SISO systems is presented. The complexity of the robust stability criterion for MIMO systems makes extension of the direct test of the elliptical bounds to M I M O systems difficult, and such extension is left for future work. The second stage in the proposed methodology is to use the identified model and uncertainty description to design future experiments. The focus here is to design control-relevant experiments, i.e. experiments that generate data most useful for developing a high-performance model-based controller. We adopt a statistical framework in which future experiments are designed based on past information. The test signal is set equal to the argument that minimizes the conditional expectation of a closed-loop performance index. This paper is arranged as follows. In the following section, the problem of deriving frequency domain bounds is addressed. 'Soft' bounds are derived in a manner similar to that of previous studies 12. In 'Robust stability tests', the problem of testing robust stability
•
Estimating uncertainty. Since model error causes controller performance to deteriorate, it is important to have a useful uncertainty description for predicting stability and performance. For example, frequency domain bounds on the system frequency response would be useful in /z-synthesis and /zanalysis. • Shaping uncertainty. It is intuitive that the test signal used in identification affects the model error distribution and, hence, the controller performance. Ideally, the test signal should be used to influence the error distribution in an optimal (for control) way. • Accounting for uncertainty. Once an uncertainty description is provided, it should be used to analyse controller performance and possibly to improve controller performance. The problem of integrated identification and control covers a broad range of disciplines, and by no means does this paper represent a complete solution. In this paper we will address the two most important stages in an integrated identification and control methodology: *To w h o m c o r r e s p o n d e n c e s h o u l d be a d d r e s s e d . Fax: (205)844-2063; e-mail: jhl!a e n g . a u b u r n . e d u
431
Integrated identification and robust control." B. L. Cooley, J. 14. Lee
432
__~
,,st
Uncertainty Quantifier
Identification Algorithm
YES P
Analysis
NO Control-Relevant Test Input Generator Figure 1 Integrated identification and control methodology
using the derived bounds is approached. It is shown that the derived bounds can be directly tested for singleinput, single-output (SIS®) systems, but only conservative tests are available for multiple-input, multipleoutput (MIMO) systems. In 'Performance tests', the problem of estimating performance based on parameter estimates is addressed. In 'Experiment design', a method for using prior information to design the test signal in a control-relevant manner is developed. Finally, 'Numerical example' presents an example illustrating how the pieces of the proposed integrated identification and control methodology fit together.
Estimating uncertainty Deriving time-domain bounds The underlying system is assumed to be described by a finite expansion in a set of orthonormal basis functions
fk(q):
fk(z) = ~ z
(1 - az) k ~ f a ~ i [a[ < 1
(5)
Other basis sets exist, but the two mentioned above are perhaps the most widely used. Note that model (1) makes sense only for stable systems. Remark: The assumptions on v(t) are not as restrictive as they may seem. If v(t) can be described as v(t) = H(q)e(t), where e(t) is a normally distributed i.i.d, random sequence of zero-mean and covariance I,y with H(q) known or estimated, the model can be put into the form of (1) by filtering the output with H-~(q). This is known as 'prewhitening'. The estimated model then has the form H-l(q)G(q), where G(q) is the system transfer function. Given (1), an estimate of ® can be obtained via least squares from a set of N input-output data. N
(9 = P Z
¢(t)yr(t)
(6)
t=l n
y(t) = ZOgfk(q)u(t) + v(t) = ore(t) + v(t)
(1)
-1
k=l
(7)
®r = [01 ... 0~]
(2)
q~(t) = [fi (q)ur(t)...f~ (q)ur(t)] 7"
(3)
In the above, u E R n" is the system input vector, y E R'y is the system output vector, Ok E R ny×n" are the (realvalued) expansion coefficients, and v C R "~ represents the effects of noise and disturbances. The sequence {v(t)} is assumed to be a Gaussian independently and identically distributed (i.i.d.) random sequence of zeromean and identity covariance. The bases fk(q) in the forward shift operator q (i.e. q-lu(t) = u(t - 1)) have the following orthonormality property:
fk(d'°)fi(e-J°~)dw =
0 k ¢ l
(4)
--Tr
Examples of such bases include fk(z)=z -l impulse response models) and Laguerre bases:
(finite
If v(t) ts uncorrelated with ~(t), then O will be an unbiased estimate of ® (i.e. E® = ®). Furthermore, t~ has the following covariance: ng blocks
Is) O" In the above, vec(.) is an operator that stacks the columns of a matrix in a vector:
vecIa'a211=[a'
Integrated identification and robust control.'B. L. Cooley, J. H. Lee Deriving time-domain bounds
Then for w = Ax C R"", nw<_nx, A c R ""×'' with full rank,
'Soft' uncertainty bounds may be derived using the X2 distribution: bt = vec(O - o)rE(elvec((2) -- O) --~ X2Co)
(10)
where p = n ny n, is the number of parameters. Equation (10) simply states that bt has X2 distribution with p degrees of freedom. This fact enables the construction of confidence bounds:
P[bt < X2×(p)] = y
433
wT(AZAT)-Iw
T T T gkl(ATSkiZoSklAT)-lgkl < X2y(p)
= FRe(dk'-Gk') ] gkl
Deriving frequeno,-domain bounds
(16)
1
Proof" See Wahlberg and Ljung 4. Using Lemma 1 with ( l l ) and (14), frequency response bounds can be derived as follows:
(11)
In the above, 7 is the probability (or confidence) that the true parameters lie within the specified ellipsoidal set, and X2(p), which depends on 7, determines the size of the set. For example, if 7 = 0.95, then X~(2)= 5.99 (see Papoulis3). Equation (11) provides a simple method for bounding the time domain parameters.
<
LIm(Gkt
(17)
(18)
Gkt) J
The above inequality describes confidence ellipses in the Nyquist plane. The ellipses are centered at Gkt and have axes defined by A TSkl Z (0Sklr A rT"
Remarks." 1. The bound (17) is based on the joint distribution of the model parameters. In other work 12, bounds are drawn based on the marginal distribution of
gr:
To test the stability of a particular controller, it is convenient to translate the time domain bounds into the frequency domain. Notice the following
Gkt(e/', (0) = A(d°')S~tvec(O)
(121)
where Gkl(d% O) is the k,l element of the frequency response matrix, and Skz is a selector matrix that selects the appropriate elements of 0 to construct the k,l element. A,(e #°) has the following form:
A(d '°) ~ [A, (d'°)... A.(d°')]
(13)
where Ak are the frequency responses of the corresponding orthonormal basis function. For example, if fk(q) = q k then Ak(# '°) = e -jk'°. Notation: In the following the arguments #~o and 0 will be dropped when it will not cause confusion. The frequency response of the underlying system will be referred to as G = G~d', ®), and that of the estimated model as d = G(d °', 0). It is useful to split the frequency response into real and imaginary parts:
r-e/<,l = [ Im(Gkt) J
Im(A)
Sklvec(O)
(14)
AT
1" A T T gkl( TSklEOSktAr)
-1
~ gkt <_ X•(2)
(19)
The result (17) suggests that the above bound actually has confidence V' such that X~/(P) = X2(2). If 7 is significantly less than 1 (e.g. V= 0.95) or n is large, the set defined by (19) can be considerably smaller than required by (17). 2. Using the bounds in (17) causes the loss of information about the correlation with respect to both frequency and the different elements of G. In the context of a robust stability test, losing the frequency correlation introduces no conservatism since the test need fail at only one frequency. However, losing the correlation between elements can lead to very conservative results since the multivariable Nyquist criterion will strongly depend on the structure of the uncertainty between different elements of the frequency response matrix. While robust stability of SISO systems can be tested using the bounds in (17), robust stability of MIMO systems depends on all the elements at a particular frequency. Using (12) and Lemma 1, it is straightforward to show:
gT(sN(.)ST)-'g < X2×(p)
(20)
where where Re(.) extracts the real part of a complex number and Im(.) extracts the imaginary part. To translate the time domain bounds to frequency domain bounds, the following result 4 is needed: Lemma 1: Let x E R "~, Z > 0 C R n*xn." and xrE Ix< 1
(151)
I Re(Gll - G l l ) Im(Gll. - Gll) g
Re(A)&l Im(A)Sli
(21)
S = [ I m ( G .... -- a < . . )
Re(A)Sn,n. Im(A)Sn,..
434
Integrated identification and robust control. B. L. Cooley, J. H. Lee
Element-by-element circular bounding.
The bounds in (20) (or even (17)) are difficult to use for M I M O systems. For use with standard analysis and synthesis tools, the bounds must be approximated with circles. It is instructive to examine the singular value decomposition (SVD) of A Skiv A r. r"
As an alternative, it is possible to retain information about the structure of the uncertainty in the frequency response elements for a given output by noticing that: AGrowk = ( ( G -
TSkl]~0
TTT : !U1.~ f ffl,k'o f2,k'O] u~TI ATSkI~aoS~IA Ukl
~
(22)
J
(23)
Idk, - Gktl <_ Cx~Co)m,~, Define
I Ac'kl 1 " A
1 E~/2(AGrowk)T
(24)
L &.,k..J
In the above, (A)rowk is the kth row of A (see (25)), and ® represents the Kronecker product• Note that A (not At) appears in (27). The following bound may be derived:
(28) In the above, AG is the complex conjugate of AG. This bound, although it incorporates information about the uncertainty structure, is also conservative. This is easy to see since, for a SISO system, the bound approximates an ellipse (17) with a disk ((28) is a magnitude bound on a complex scalar for SISO systems). Unfortunately, the disk is larger than the bound defined by (23). However, (28) contains information about the relationship between the frequency response elements for a given output. Consider the following singular value decomposition: (a ®
AGrowk ~ [Gkl -- Gk,... dk,~. - Gkn,]
O'l,kl
0
-"
.
"
(27)
AGrowk((A ®ln.)P(A* ® I,~.)) l(AGrowk)T < X2v(p)
Gkl
where al,kl ~ Cr2,kl. Thus, lies in an ellipse of center (~kl with axes ul and u2, of lengths ~ and ~v/~,~/. The smallest circle that overbounds this set has radius (see Figure 2). This circle is described by
Ac,k =
G)T)rowk(AT @In.)
(25)
0
I,,.)P(A* ® i,,.)
(26)
(29)
where * denotes the complex conjugate transpose, Us is unitary, and E~ is diagonal. Define
As,k =
~ac, k
U, sG
=
I As'kl " ]
A
L A,'k.°
l
~l/2G(AGrowk) T
(30)
"
0
0
-..
0
(71,kn.
It is easy to show via (23) that
]Ac,kzl
< 1, l = 1. . . . . n,.
Row-by-row structured bounds. Obviously, approximating the bounds of (17) with (23) is conservative when ~rl,kl > ~rZ,kt since the circular bounds include values for Gk/ not in the elliptical bound• However, a greater criticism of using these bounds is that they are independent bounds on each element of the frequency response matrix and therefore do not incorporate information about the structure of the uncertainty. This is potentially very conservative. Im
In the above, As,kt, I = 1. . . . . n, are complex random variables. It is easy to show via Lemma 1 that (28) implies that IAs,kll _< 1, l = 1. . . . . n, Remarks. • Both (23) and (28) are conservative bounds. However, the relative conservativeness of these bounds depends on the characteristics of Eo. If this matrix is such that there is little correlation between elements of the frequency response matrix, then it is expected that (23) will be less conservative than (28). On the other hand, when such correlation exists, (23) can be very conservative in comparison to (28). In general, both bounds should be computed, and the structured singular value should be taken to be the minimum value for the two bounds.
O-lU 1
Robust stability tests SISO systems
Re Figure 2 Overbound approximation to (17)
While it would be simple to use /z-theory (i.e. existing tools) by approximating the elliptical bounds in (17) with over-bounding circles, such an approach would clearly lead to conservative results. Instead, we will
Integrated identification and robust control." B. L. Cooley, J. 14. Lee derive a robust stability test that explicitly uses the elliptical bounds by translating the stability condition for SISO systems into a test based on the elliptical bounds. For SISO systems, the robust stability criterion simply asks whether G C = - I for any model in the model set and a given controller C. This simple criterion is easy to test upon considering the following representation of the product of GC: G C -- Re(C)Re(G) - Im(C)Im(G)
(31) + (Im(C)Re(G) + R e ( C ) I m ( G ) ) j
!Re(c) Im(C)
1 rRe(oc)l LIm(GC)]
Im(G)
Re(C)
=
(32)
Note that the (frequency-dependent) transformation matrix CT maps from the complex plane to the complex plane. Define
g'
A [Re((G-G)C)] O)C)
(33)
435
ellipsoidal set centered at G C with axes determined by (CTATE(.)ATCTT) that does not contain G C = - I is determined by T
gr
T T (CTAT]~(')ATCT)
I
gc
<
rain h
(39)
Remarks. •
1. Determining the maximum confidence is very desirable since Y is usually chosen somewhat arbitrarily. Knowledge of the maximum Y would allow the engineer to make an intelligent choice between accepting the model or collecting more data to refine the model. For example, if the acceptable confidence is 7=0.95, then the engineer might decide not to retest if it were determined that min~ b = X~.94(F/). Ideally, the y associated with minoj b could be monitored on-line and used to decide when to stop testing. 2. The robust stability test outlined above is similar in concept to that of a number of previous papers 5 7. Note that the above robust stability test and the tests described in the cited works are exact tests (i.e. no approximation of the bounds is required). M I M O systems
The mapped bounds (n,, = ny = l) become gc Zx, g' <- x?(n)
A
~g, =
(34)
T
T
(35)
CTA T~(.)A TCT
To check the bounds for robust stability, define = L
Im((~C)
z_~ g'
Re(dC) + l Im((~C)
(36)
Since and (~ and C are known, b is a known function. Robust stability can be determined by checking the following inequality at each frequency: 9
b < X~,(n)
(37)
If (37) holds at any frequency, then the uncertainty set defined by (34) contains the point G C = 1 in its Nyquist band, and robust stability does not hold with confidence % The above stability test can be extended to determine the maximum confidence for robust stability. Notice that b >_ m i n h
(38)
2 Thus, if y is selected such that X× < min~ b, then (37) will fail at every frequency. In addition, the largest
To avoid conservativeness, MIMO systems should be tested in a manner similar to that outlined in the previous section for SISO systems. Unfortunately, the multivariable Nyquist criterion, which asks whether the map of det(& + GC) contains 0 for any model within the model set characterized by the bounds in (20), is complex, and it is difficult to extend the simple technique developed for SISO systems to MIMO systems. As an alternative, structured singular value (la) theory can be used to test stability./~-Theory is a powerful tool for checking robust stability, but the test assumes that the frequency domain bounds are disk-shaped. While the bounds of (17) cannot be used directly, those of (23) or (28) can be used. Since the relative conservativeness of the two bounds depends on the structure of Z(.}, a stability test should be performed for each type o[" bound, and the least conservative result (smallest values of the structured singular value) should be used. The bounds in (23) can be used as lbllows. First, uncertainty descriptions can be constructed for each output. It is apparent from (24) that
AGro,~k -- [1-
0
.11][{£ ." ,6,j,i
•.
0 ]
•
0
0
A,. k.,,
(4o)
Integrated identification and robust controL B. L. Cooley, J. H. Lee
436
The uncertainty descriptions for each output can be combined as shown in Figure 3 with Z = Zc (see below) and D = I,,,y. This system can easily be put in standard M - A form 8 with
M = Wl,sC(Iny -~- G C ) - I W2
,,,',,
M = Wl,cC(In~ +
2
(~C)-I W2
"½
f/-}
=
,,, blocks
with the additional definitions
ny blocks; ",
Us zx
".°
"..
"•
"..
(42)
il
"..
• "•
1!;
•,
"•
•
"..
". •
0
...
0
A
0z~!0
n. elements
0
''"
,.
*,
••"
0
...
0!
(43)
0
(44) Ec.,~
For the bounds of (28), the uncertainty description for a particular output Yk can be written as follows: Akl
0
...
0 ]
..
".
".
".
0
0
Akm
Aarowk = [ 1 . . - 1 ] 0
x
•
" "
1
(45)
V/ 2(p)Zs1/2 UsT
The uncertainty each output and = Zs and U = in standard M - A
descriptions can be constructed for combined as shown in Figure 3 with Us (see below). This can easily be put form with
.
The uncertainty description for the M I M O system
0 ]
~
1
(49)
A structured singular value test can be conducted for M of the form in (41) and for M of the form (46). Since both tests are based on conservative bounds, the smallest structured singular value obtained from the two tests can be used as an upperbound on the structured singular value that would result from testing the true ellipsoidal bounds• R e m a r k : In 'Numerical example' an example will be presented to demonstrate the potential conservatism of the robust stability tests based on (23)/(41) and (27)/ (46). It is important to note that such conservatism arises from force-fitting the ellipsoidal bounds into bounds that can be used in standard/z-analysis. There are two ways to reduce (or even eliminate) the conservatism• First, less conservative approximations of the ellipsoidal bounds may be sought. Second, alternative robustness tests that directly use the ellipsoidal bounds may be developed. It is our belief that the latter approach is the correct one, and it is the approach we used for SISO systems (see 'SISO systems')• However, for SISO systems, the stability test involves only checking whether a point lies within an ellipse• For MIMO systems, the stability criterion is much more complicated (det(I,~ + G C ) = 0) and not so easily checked. Clearly, the search for such alternative robustness tests deserves a paper in its own right and is beyond the scope of the current paper.
Performance
Figure 3
...
n~ blocks
n. elements
0
(48)
0
0
$ ----
y
]~e,l
.•
Z~
with the additional definitions ...
i1 .•
ny b]Cocks
ny blocks
iS[1
(47)
(41)
I"]} •
(46)
tests
While robust stability is desirable, the ultimate goal is good performance. For example one might be interested in obtaining good mean-squared error performance: J = EeT(t)e(t) = tr(E{e(t)er(t)})
(50)
Integrated identification and robust controk B. L. Cooley, J. H. Lee In the above, e represents the difference between the reference signal and the system output (see Figure 4). It is straightforward to show
437
An interesting result is that minimization of the righthand side of (54) is equivalent to minimizing the weighted trace of the parameter covariance matrix: min J ~ min W'Zo
e=S(Iny + (G - (~)(}-I/})-1 (r - d)
(56)
(51) A formula for the weight W can be found elsewhere 9.
Remarks." In the above S = (I.y + GC) -I and H = I . , - S. d is assumed invertible for convenience. Using the above expression and Parseval's relationship, it is easy to show that (50) becomes J = t r ( I~
f~ E{S(I< +(G-O)G
1/})-1 (52)
1. The controller does not explicitly aflpear in (54). Instead, the sensitivity function S is specified. Since the model G is available, determining S for a given C is straightforward= Note that it is a common practice to specify S and translate it to C once G is obtained. 2. Care should be taken in using (54) since it only holds when the closed-loop system is stable.
x q%(/< + (G - G)d-l/-})-*S*}dw)
Experiment design In the above, the superscript * represents complex conjugate transpose (A-* is the inverse of A) and qbn is the spectrum of rl. Equation (52) directly measures the expected performance of the controller. However, it is difficult to evaluate because the unknown plant G enters into the expectation in a nonlinear fashion. When the plant estimate is close to the true plant, (51) may be linearized about (~, which is the best estimate of the plant: e
+
(53)
-
Thus, for unbiased (~
(54)
The stability test and the performance test described in the preceding sections provide criteria for accepting or rejecting a model. When the model is rejected, more data must be collected to refine the model. The information gained from previous experiments can be used to extract more control-relevant information from the system through test signal design. It is intuitive that the optimal experiment design will depend upon the system dynamics. Thus, it is necessary to employ an iterative scheme where the design is computed based on a prior distribution of the plant parameters. The test signal design problem is formulated as follows. First, assume that a prior distribution on the plant is given in the form of a nominal model (2)k, and its covariance E(.)k. This prior might be specified by the engineer, or it might come from the estimates provided from a previous experiment. When new data are collected, the model can be updated as follows:
N
) 1 (57)
The second term in the integrand of (54) is a linear combination of E{((~ - G)*(G - G)} which is given by Ok+l = Z (Ok -1 + I E{(G - G)*(G - G)} =
(P(t)yr(t)+ E6~It2)k
(581
ny(A r ® I..)*P(A r ® I..)
(55) Notation: Here, the notation is extended to incorporate In the above ® represents the Kronecker product.
r
_
~
the iterative design by defining, for example, G~-G(d°~,~)k). This is not to be confused with GrowkThe test signal is designed to minimize the conditional expectation of the performance index in (52) subject to any constraint on the inputs:
y u~kP+tl= arg min
Figure 4 The closed-loop system resulting from identification of model G and constructionof model-basedcontrollerC
Uk+l E U
E{JlOk, I2~.~}
(59)
(60)
Integrated identification and robust controk B. L. Cooley, J. H. Lee
438
In the above, Ug+l is the test signal used in the (k+ 1)st experiment. U is a set that expresses constraints on the test signal and is usually characterized by magnitude and rate bounds. Note that J depends on the plant, the reference signal/disturbance spectrum, and the model obtained after the experiment: J
=
tr
Sk+l(/., + (G
x O,(I,> + (G
-
X S~+l](~kZOk
-
--
^
^-I ^
Gk+l)Gk+lHk+l
)-I
Gk+l)Gkll/~k+l) -*
})
doo
(61) The plant and the model are unknown and treated as random variables. The minimization (59) is difficult to solve because the objective function contains the unknown plant and the model, which enter in a nonlinear fashion. It is convenient to linearize the error expression, this time with respect to both the plant and the model. The linearization is taken around the best estimate of the plant, Gk: e ~ (Sk+l + Sk+l((}k+l - G)G/1/Qk+I)~
(62)
Note that, in general, Sk+l and /}k+l depend on (~k+l. However, if a direct synthesis approach to designing the controller is taken (i.e. the sensitivity function is specified), then Sk+l and Hk+l become design variables (i.e. their values are specified by the engineer). The minimization in (59) becomes U~kpt' = arg min tr
~
E Sk+l (dk+l - G) -Tr
^
1
^
^
*
^
x G k g k + l ~ ( G k I-I~k+l
(63)
-- Gk+l) S~+I]OkZO k d o
It is straightforward to show that the above objective function is simply a linear combination of the elements of E(() - G)*(G - G), and can thus be related to the test signal through (55). As a result, the above optimization can be solved, but the potentially large dimension of Uk makes computation an issue. A suboptimal way to compute this optimization is to divide the test signal into small portions and solve for each portion. Furthermore, (63) is equivalent to minimizing the weighted trace of the parameter covariance matrix: min J ~ min WZo
resulting model (and the uncertainty) by denoting it with an index corresponding to the experiment being performed (i.e. &+l). Inclusion of uncertainty in the controller design greatly complicates the experiment design, however. Note that, in this case, Sk+l is not known prior to the experiment and hence is a random quantity. To solve the problem rigorously, one must embed into the expectation the dependence of S upon the future model and the uncertainty. This topic is worthy of a paper in its own right, and we defer an indepth study to future research. In practice, since the objective is only^to improve an experiment design, assuming Sk+l = Sk will generally be sufficient.
(64)
A formula for the weight W can be found elsewhere 9. Remark." One could potentially consider a more complicated robust control system design method, rather than a direct synthesis method (for which the closedloop functions are fixed). We implied that the desirable closed-loop nominal response may depend upon the
Numerical example Preliminaries In this section, we will test the integrated identification and control methodology discussed in this paper on a nonlinear high purity binary distillation column model. The model for the column has been set forth previously l° and has been considered in a number of recent works 11,j2. Briefly described, the column consists of 25 trays, a reboiler, and a total condenser and is designed to separate a binary mixture of methanol and ethanol into product streams of 99% purity. The two controlled variables are the temperatures on trays 21 (T21) and 7 (TT). The two manipulated variables are the reflux flow (L) and vapor boilup (V). The nominal values are T21=614.1°R, T7=636.8°R, L-3.3835mol/min, and V = 3.856 mol/min. The column model includes tray-bytray energy balances and tray hydraulics. The objective in this example is to use the integrated identification and control methodology to identify models for use in building controllers. The controllers are specified by designing a desired closed-loop response. Here, the complementary sensitivity function is fixed so that the nominal response of the outputs to the reference signal is first-order with time constant 4: /} = (1 - e-1/4)z -1 1 e_l/4z_ 1 12
(65)
This corresponds to the following controller: C(z) - (1 - e-1/4)z -' ~ 1 1 -
z -1
(z)
(66)
The above controller is used when the model contains no zeros outside the unit disk. When the model contains zeros outside the unit disk, the model is separated into all-pass and a minimum phase portions as follows: G(z) = GA (Z)GM(Z)
(67)
where Ga(z) is all-pass, and GM(Z) is minimum phase. See Morari and Zafiriou 8 for details on this inner-outer factorization. The following controller is used:
Integrated identification and robust control." B. L. Cooley, J. H. Lee
C(z) = (1
-
15
e 1/4)2-1
i --- ~
439
GMi (z)
'l
(681) 10
While the above controller may give a different nominal sensitivity function from the one used in the experiment design, the sensitivity functions are usually simillar. The above controller will maintain closed-loop stability when the model contains unstable zeros, however. Note that (66) and (68) are controllers with integral action. To obtain the models, preliminary experiments are performed to determine the settling time (around 240 min). Based on these results, we decided to sample the column every 4min and model it using a finite impulse response model consisting of 60 FIR coefficient matrices. Two sets of experiments are performed on the column. In the first set, 1000 samples are generated using an independent random binary signal (RBS) for each input channel of deviation +0.012mol/min from the nominal values. In the second set of experiments, the experiment design methodology of 'Experimental design' is used. The experiments are designed by breaking each experiment into four sections, each consisting of 250 samples. The first 250 samples are generated using an independent RBS for each input of deviation +0.012mol/min from the nominal values. A preliminary model is then identified based on this data, and this model is used to design the input signal for the next 250 samples with each input constrained to a maximum deviation of +0.012mol/min. The model is then updated recursively, and the experiment design procedure is repeated using the new model. A total of 1000 samples are collected (250 RBS, 750 designed). For each experiment, the outputs are corrupted with a zero-mean independently and identically distributed random sequence of covariance (0.1°R)212. Controllers are designed based on the resulting models as described above.
Analysis of stabilio, Uncertainty bounds can be constructed for the identified models as described in "Estimating uncertainty' and used to analyse the stability of the closed-loop system for a given controller. Plots of the structured singular value (/x) for the two uncertainty tests are given in Figure 5. It is interesting to note that (23) provides a tighter /~ than (28) for the RBS trials, while (28) provides a tighter/x than (23) for the optimal design. This suggests that it is crucial to capture the structure of the uncertainty when computing the stability measure. One particularly notable aspect of this trial is that the /x-analysis result seems to indicate that the closed-loop system should frequently be unstable, both for the RBS experiments and for the optimal designs. Judging from Tabh, 1, this may be true for the RBS designs, but it is not true for the optimal designs. Moreover, it is difficult to discern from the figure which design is better, despite the fact that the simulation results proved the optimal design to result in a stabilizing controller more often.
iI
1
t~\//"A'Xv/'41, 0
i'L
fi ~ ,~, , ,A'
~' ~
~,
I l, •
0!f
012
013
014 frequency
0'5
016
0L7
, 01
012
0~3
0i4 frequency
015
i 06
i 07
25
>15
(b)
10 5
O0 25 Ii
2o 1: > 15
!o
L,:
/t
,'
" "
1
5 ji.L/~tl_t'~ll ,'/i~izt/~ x lit k' 11 ,1' , 01
02
03
,i /
04 frequency
,.
05
!'
^ 06
11~ , ,I 07
Figure 5 typical result: (a) structured singular value (ix) for RBS design (solid line, computed from (23)/(4]): dashed line, computed from (2?)/(46)); (b) ix for optimal design (solid line, computed from (41); dashed ]me, computed from (46)); (c) comparison of the RBS and optimal design (solid line, smallest IX for RBS design: dashed line, smallest ls for optimal design)
This indicates the conservativeness of the robust stability tests described in 'MIMO systems'. Clearly, these results motivate the search for an alternative robust stability test that can use the derived ellipsoidal bounds without approximation. One approach for reducing the conservatism in MIMO robustness tests is that of highly structured stability margins, developed previously 13"14 While the concept is applicable to the problem developed in this work, algorithms that explicitly deal with ellipsoidal bounds are not developed in these studies. The development of such algorithms is beyond the scope of this paper and is left lk~r future work.
Comparison q/ RBS and optimal design It is of interest to compare the performance of controllers designed based on the models identified from each experiment. Table 1 shows the results of the two sets of experiments. Clearly, the experiment design technique greatly improved the closed-loop behavior, both by yielding more stabilizing controllers and better performing controllers. Note that the unstable runs are such that the matrix G(0)(~ 1(0) has eigenvalues with negative real part (G(0)) is the steady state gain matrix Table I
Simulation results
No. unstable responses Average ISE (stable)
RBS design
Oplimal design
45 {).7827
8 0.3858
Integrated identification and robust control." B. L. Cooley, J. H. Lee
440 614~
and performance. If it is determined that the model must be refined, the information gained from previous experiments can be used to help design future experiments so that more control-relevant information is collected.
6143E 614
~ (a) ~
614 25 814.2
~'614 15 614 1
References
614,05 614 0
I 100
J 200
I 300
I 400 time
I 500
L 300
L 700
J 800
500
600
700
800
837¸1
307 '"
"~"~~
-
~
~ 3 0 6 g5 {b) r,. 636g >, 636 85 636 8 638 75
0
1 O0
200
300
400 IJrne
Figure 6 (a) T21 output response (solid line, RBS design; dashed line, optimal design; dotted line, setpoint); (b) T7 output response (solid line, RBS design; dashed line, optimal design; dotted line, setpoint)
based on the linearized model). In such cases there is no controller (based on the model) with integral action that can stabilize the closed-loop system 8. The results from a particular trial are shown in Figure 6.
Conclusions This paper develops an integrated identification and control scheme that provides not only frequency domain uncertainty bounds but also tests for stability
1. Goberdhansingh, E., Wang, L. and Cluett, W. R., Robust frequency domain identification. Chem. Eng. Sci., 1992, 47(8), 19891999. 2. Goodwin, G. C., Gevers, M. and Ninness, B., Quantifying the error in estimated transfer functions with application to model order selection. IEEE Trans. Automatic Control, 1992, 37(7), 913928. 3. Papoulis, A., Probabilty, Random Variables, and Stochastic Processes. McGraw-Hill, New York, 1991. 4. Wahlberg, B. and Ljung, L., Hard frequency domain model error bounds from least-squares like identification techniques. IEEE Trans. Automatic Control, 1992, 37(7), 90~912. 5. Tsypkin, Y.-Z. and Polyak, B. T., Frequency domain criteria for LP-robust stability of continuous linear systems. IEEE Trans. Automatic Control, 1991, 36, 1464-1469. 6. Blernacki, R. M., Hwang, H. and Bliattacharya, S. P., Robust stability with structured real parameter perturbations. IEEE Trans. Automatic Control, 1987, 32(6), 495 506. 7. Guzzella, L., Crisalle, O. D., Kraus, F. J. and Bonvin, D., Necessary and sufficient conditions for the robust stabilizing control of linear plants with ellipsoidal parametric uncertainties. In Proc. 30th CDC, Brighton, UK, 1991, pp. 2948-2953. 8. Morari, M. and Zafiriou, E., Robust Process Control. PrenticeHall, Englewood Cliffs, NJ, 1989. 9. Cooley, B. L. and Lee, J. H., Control-relevant experiment design for multivariable systems. Automatica, submitted. 10. Weischedel, K. and McAvoy, T., Feasibility of decoupling in conventionally controlled distillation columns. Ind. Eng. Chem. Fundamen., 1980, 19, 379-384. 11. Li, W. and Lee, J. H., Control relevant identification of ill-conditioned processes. Comput. Chem. Eng., 1996, 20, 1023-1042. 12. Chien, I.-L. and Ogunnalke, B., Modeling and control of highpurity distillation columns. In AIChE Annual Meeting, paper 2a. Miami, 1992. 13. Figueroa, J. L., Desages, A. C., Romagnoli, J. A. and Palazoglu, A., Highly structured stability margins for process control systems. Comput. Chem. Eng., 1991, 15(7), 493-502. 14. Rotstein, H., Desages, A. and Romagnoli, J., Calculation of highly structured stability margins. Int. J. Control, 1989, 49(3), 1079-1092.