A Frequency Domain Criterion for Feasible Parameter Set Recursive Bounding

Copyright (i:) IFAC System Identification, Copenhagen, Denmark, 1994

A FREQUENCY DOMAIN CRITERION FOR FEASIBLE PARAMETER SET RECURSIVE BOUNDING * G. Zappa o ,

A. Vicino"

and

A. Tesi o

o Dipartimento di Si$temi e InJormatica, Univer$itd di Firenze Via di S. Marta 3, 50139 Firenze, Italy Facoltd di Ingegneria, Univenitd di Siena Via Roma 56, 53100 Siena, Italy

00

Abstract. This paper presents a new procedure for outbounding recursively the feasible parameter set of a linear model through parallelotopes. The key feature of the algorithm is the uptdating rule based on a frequency domain optimality criterion. The parameter set estimate is updated by requiring that the feasible frequency response set in the complex plane is "minimal" in some sense. The updating law is designed so that among all the estimates with "minimal" frequency response set, the minimal volume estimate in the parameter space is selected.

Key words. Recursive identification; parallelotope; parametric uncertainty; set membership identification; linear model.

1. INTRODUCTION

sponse is negligible. This hypothesis, which is often met in real applications, allows one to split the objective of minimizing the frequency uncertainty bounds of the identified model, in two separated tasks. The first one consists in constructing an estimate of the parametric model uncertainty set ensuring that the corresponding set of frequency responses is as "small" as possible. The second objective is to find the minimum H oo bound compatible with the a priori information and the measured input-output data. In the present paper, we focus on the former objective. The main objective is to build a procedure which allows one to outbound adaptively the parameter uncertainty set through parallelotopes. The main feature of the algorithm proposed here, which differentiates it from that in Vicino and Zappa (1993), is a new updating criterion, which consists in minimizing a norm bound on the frequency response value set at a given value of the frequency.

The objective of most contributions in the literature on H 00 identification consists in deriving a nominal model for the plant and guaranteed frequency domain uncertainty bounds to be used for robust control design purposes. Most of the recent contributions refer to nonparametric approaches (see, e.g., Dahleh and Smith (1993), Milanese and Vicino (1993), Helmicki et al. (1991), Gu and Khargonekar (1992), Chen et al. (1992», while few are oriented to mixed models where uncertainty is described through structured parametric perturbations and unstructured H oo uncertainty. Moreover, except for few cases (Younce and Rohrs, 1990; Kosut et al., 1992; Wahlberg and Ljung, 1992; Vicino and Zappa, 1993), almost the totality of the existing contributions provide one shot (off line) solutions, which are not very suitable for robust adaptive control schemes. In the present paper we take a recursive approach to the robust identification problem in the same line as Younce and Rohrs (1990), Kosut et al. (1992), Vicino and Zappa (1993). More specifically, we assume that the system model accounts for both parametric perturbations and unmodelled dynamics. The basic underlying assumption is that the parametric model structure gives a reliable representation of the system in the low frequency range, while the unmodelled dynamics uncertainty is concentrated in the high frequency range, where the parameteric model frequency reo

The paper is organized as follows. Section 2 introduces notation and some preliminary results used in the following development. Section 3 provides the main results of the paper on the approximation of the feasible parameter set. Finally, some comments and concluding remarks are reported in Section 4.

2. NOTATION, PRELIMINAR RESULTS AND PROBLEM FORMULATION Suppose that a set of data is available (u(k), y(k»,

Work supported by Italian MURST and CNR.

5J5

k = 1, ... , N ,

where IW(ej"')1 ~ O. The second objective is to estimate the minimum H oo bound (6, quantifying the magnitude of the nonparametric perturbations in the high frequency range (frequencies at which IW(d"')1 ~ 1) .

where u and y represent the measured input and output of the system to be identified. It is assumed that the data is produced by an unknown linear, time-invariant system Go(q) (q represents the unit forward shift operator), and that the measured output data is corrupted by a noise term e, I.e. ,

y

Let Y(N) = [y(I), ... , y(N»)' and let ~(N) denote an (N, n) matrix depending on the input vector U(N) [u(I), ... , u(N»)' according to the specific linear parameterization chosen (for example, if the model is simply a FIR filter, ~(N) is a lower triangle Toeplitz matrix). Denote by
= Go(q)u + e ,

=

where e is bounded according to an loo norm, i.e.,

For our identification problem, we will assume that the true system admits the following representation

Go(q)

= G(q; 0) + 6G(q)

(2)

,

=

where (l:, k 1, ... , N depend on the measurement error bound (. From (2) it turns out that eN is a convex polytope.

where 0 E Rn is the vector of unknown parameters, G(q; 0) represents the parametric component of the model, while 6G(q) stands for possible nonparametric uncertainty accounting for un modelled or difficult to model dynamics. Several different hypotheses can be done on G(q; 0) and 6G(q). Regarding G( q; 0), we will assume a linear parameterization, i.e.,

G(q;O) = O'g(q) ,

Let us now introduce some terminology on polytopes. Denote by V E Rn a bounded convex polytope, i.e., the solution set of a system of linear inequalities

V == {U ERn: AO ~ I} ,

(1)

being A E Rrxn full column-rank. The hyperplane

where g(q) = [gl(q), ... ,gn(q»)', being gi(q), i = 1, ... ,n elements of an assigned basis in the model space (e.g., Laguerre or Kautz filters, or gi(q) q-i if a FIR model is assumed). As far as 6G(q) is concerned, we assume that

=

6G(q)

(3)

(1

== {O E Rn : pO = I},

p' E Rn

(4)

is called a supporting or tangent hyperplane for V if it has non void intersection with V and

= W(q)~(q) ,

pO

~

1 or pO

~

1,

'rIB E V.

(5)

The set

where W(q) is a rational proper, stable weighting function and ~(q) is an H oo norm bounded perturbation

S(p,c) == {O: IpO -

cl

~

I}

(6)

is called a strip.

1I~(q)lIoo ~ (6 .

Definition 1. - The strip S(p, c) is said to be It is worth while to make few observations on the role played by the weight W(q) in the overall frequency domain oriented identification approach considered in the paper. The function W(q) is assumed to be known. Although its shape depends on the specific case at hand, W(q) will be in general a high pass filter. The main objective for introducing this type of filter is to separate the frequency range where the parametric model provides a reliable representation of the system from the range of frequencies where all information on the system is lumped in the H 00 norm bound (6. This allows the designer to decoupie the parametric uncertainty description from the nonparametric component, so that the identification objective can be split in two different goals. The first one is to derive parameter uncertainty set estimates which minimize the area of the model frequency response set at each frequency w

tight with respect to V if both hyperplanes bounding the strip are supporting hyperplanes for V . • Clearly any polytope given by (3) can be represented as the intersection of m (m ~ r) tight strips, i.e., m

V = {nS(Pi,Ci)} .

i=l For notational simplicity, we will denote S(Pi, Ci) by Si whenever no ambiguities arise. If m = n, the polytope n

is called a paraIIeIotope. In this case

V = {O: IIpO - clloo ~ 1} == P(P,c) 516

(7)

where P

==

[p~,

M(U) is a convex polytope, too. In addition, if U is a parallelotope, M(U) is a parpolygon, i.e., a

... , P~)' E R nxn e == [Cl, ... , en]' E Rn

(8) Notice that, since P(P, c) is bounded, P is invertible. It is simple to compute the vertices of P(P,e). In fact, by setting T = [tl,'" ,t n ]

== p- l ,

()e

== Te,

we get from (7) P(P, c) = {() : () = ()e

+ TO, 1181100

~ I}.

(9)

Hence, it turns out that the 2n vertices of P(P, c) can be expressed by

polygon with at most 2n vertices symmetric with respect to a center. Let us now state the problem addressed in the next section. We consider time sequences of polytopes and approximating parallelotopes denoted by V(k) and P(T(k), ()c(k», respectively. Of course, in the specific estimation setting delineated before, strips like (6) concurring to generate recursively polytopes are represented by parameter uncertainty sets determined by each measurement and the relative error bound. For example, the strip S(k), is given by

n

v/=()e+LQ/jt j ,

Q/jE{-I,I}, I=I, ... ,2 n

.

j=l

(10) We will often refer to P(P, c) as P(T, ()e), to indicate explicitly the dependence of a parallelotope on T and ()e. We will simply denote a parallelotope by P, whenever no ambiguities arise.

Remark 1. - Notice that S(p, c) = S( -P, -c) . •

Assume that, at each recursion k, a parallelotope P(k) == P(T(k),()e(k» and the strip S(k + 1) are given. Define the polytope V(k

+ 1) == P(T(k), ()e(k» n S(k + 1)

and denote by p. the minimal area parallelogram in the complex plane outbounding M(V(k + We look for a parallelotope P( k + 1) such that

1».

From (9) it follows that the volume of a parallelotope P(P, c) is given by the expression

2n

vol(P(P, c» == Il(P(P, e» = Idet (P)I Since the purpose of the following sections is to approximate eN through parallelotopes enjoying frequency domain optimality properties, we introduce the definition of (frequency dependent) "model value set" associated to an uncertainty set in Rn. Let Mw (.) : Rn -+ R 2 be a map defined for a fixed frequency w as () E Rn .

M(P(k

+ 1» = p•.

We report now some important results which will be exploited for deriving the main contribution of the paper. To this end, let us consider a polytope V defined by the intersection of n + 1 strips n

and introduce the n + 1 parallelotopes p(i) and relative matrices p(i), T(i) and ve~tors e(i), ()~i), obtained by eliminating the strip Si, i = 0, ... , n, I.e.,

n n

Since we adopt the linear parameterization (1), we have

p(i)

jW

Sj},

(14)

j=Oj~i

" ,n ]' , P (i) ..:.. - [P'o"",Pi-l,Pi+l"",P

where Mw E R 2xn is given by

jW

_ [ffil(W)] _ [Re 91(e ) ... Re 9n(e )] w ffi2(W) - Im 91(e1W ) ". Im gn(e1W) (11) Definition 2. - Let U be an uncertainty set in parameter space. For a given frequency w, we define the model value set Mw(U) as M

== {

r
It is easy to check that if U is a convex polytope,

= T(i)e(i) :

n

Po = L')'iPi ;=1

(12)

•

()~i)

Moreover, assume that the components of the vector Po with respect to the basis provided by {PI,'" ,Pn} are available, i.e.,

Mw(U) = {( E C : ( = Mw(), for some () E U} . In the sequel, we will drop the explicit dependence of Mw and MwO on w, whenever no ambiguities arise.

= (p(i»-l,

where, without loss of generality (see Remark 1), we assume that ')'i ~ O. In addition, let us set ')'0 1.

=

The following lemma relates the parallelotopes pe;), i= I, ... ,n, to ptO). 517

Lemma 2. - Assume that PI E span {ml' m2} with PI PIM. Then, M(P) is the parallelogram P(T, Bc) where

la, 1993) - Let

=

.

,. .

= [i l , Ej=2 cijij ] Bc = [i l , ... ,in]c

We shall prove th ists a correspondil is oriented and SI M(P nSo) = MP

T if j

=i

if j

#i

being cij

sign (it i j ).

n

•

VI

ppa, 1993) - Let the ht with respect to V. linimal volume outwhere

= Bc+0/1i!+012(Lcijij),

0/1,012 E {-1,1}.

j=2

Lemma 3 - Assume that both PI and P2 E span {ml,m2} with PI = PIM and P2 = P2 M . Then,

Proof If PI, P2 E span {ml' m2}, from (16) it fol0 for j ~ 3. Therefore, (17) implies lows that ij that M(P) M(Sl n S2) is the parallelogram with vertices VI, I = 1, ... ,4, given by

•

= =

an be tightened by Ions, Lemma 1 and >rocedure for the reolume of the uncerng the strip So pront. In the next sec.is procedure can be : the model value set

r Clearly, by constr So whenever So i: E span {ml' m2} that Po = PoM. F

Pc where PI and P2 P2 = P2 M , respe< C2 = C2 and Co = Then, we have th, Lemma 4. - Le' P ni=l Si is ori

=

=)

Recalling (see (14», th

Proof

= [i l ,i2 ] [~~] . Hence, we get M (P) = S(Pll Cl) n S(P2' C2)

p(j)

Bc

as with ( 18)

Now, since it turns out that (see (16»

M(V) = {n;'=, where the last e< for j ~ 3, M(P p(j) are oriented Let us now prov Since So and S2 applied. Hence, i

exploiting (18) we get

(16)

Consider now pC get

re given by

Hence, the result follows.

•

Definition 3. - The strip S(p, c) is said oriented according to the mapping M if P E span {ml' m2}' The parallelotope P = {ni= 1 Si} is said oriented if both SI and S2 are oriented . •

at M (P) is the con-

1, I},

III

M(P n So)

where

, expression for the tainty parallelotope ose, notice that the .s (compare (11»

_

where

0/1,012 E {-I, I},

AIN FEASIBLE PROXIMATION

.. ,q.p

•

Thus, M(P) is a parallelogram.

{Ii} }

III

Co = -

Proof From (16) it follows that Plij = 0, j ~ 2. Hence, the vectors i j are parallel for j ~ 2 and M(P) is the convex hull of the vectors VI, I = 1, ... ,4, given by

(15)

n

=

_

Po = -

Let us now consider the polytope V = (P n So) where the parallelotope P {ni=l S(pi, Ci)} is oriented. Moreover, let So = S(Po, co) with

=

I = 1, ... , 2n , (17)

n

Po = LliPi

.. ,in]c .

i=l

er which conditions

where we assume that li ~ 0 (see Remark 1). 518

Since S2 is orienl j ~ 3, applying. the parallelotop<

We shall prove that for any given So there exists a corresponding strip So = S(Po, co) which is oriented and such that the model value set M(P nSo) = M(P nSo). To this purpose, define

= 1'IPI +r 1'2P2 C = 1'1 Cl +r 1'2 2 + Co - rPoOc

Po Co where

(19)

1'IPI

= =

t<1)

Hence, M(P(1» is the parallelotope defined by the strips S(Po, co) and S(P2, C2) and, thus, (21) M(So n S2)' Using a similar implies M(P(1» reasoning for M(p(2» and taking into account that M(PCO» M(SI nS2), we get

=

= {n?=OS(pi,Ci)}

•

+ 1'2P2 r

We are now able to give the main result providing the solution of the following optimization problem. Given the oriented parallelotope 'P and the new strip So = S(Po, co) in the parameter space, find the oriented parallelotope p., including P n So, such that M(P*) has minimum area, and p. has minimum volume.

(20) Theorem 2.

=

=

- Let the oriented parallelotope

P = {ni=1 S;} and the new strip So be such that Po = 2::7=1 1'iPi, 1'i ~ O. Moreover, assume that all the strips are tight and define 1'0

= r.

i· arg{maxi=0,1,2hd} j* = arg{maxj=3"nhj}} .

Lemma 4. - Let V = P n So and assume that ni=1 Si is oriented. Then,

Then, the parallelotope P* is given by

Proof Recalling the definition of parallelotopes

where

=

pU) (see (14», the model value set can be written Si

as

Si· where the last equality follows by noticing that, for j ~ 3, M(pU» = M(PCO» since p(O) and pCj) are oriented. Let us now prove that M(PCI» = M(So n S2). Since So and S2 are oriented, Lemma 3 can be applied. Hence, it turns out

Consider now pCI) get

= P(r
__

Proof From Lemma 4 we know that

where Si = S(fJ;,ci), i = 0,1,2. From (20) we have

r·

.

.

r

.

1'2

.

.

is the minimum area parallelogram containing

M(V), i.e., J.l(,P*) ~ J.l(M(P»

=

Since S2 is oriented and, by definition, t)O) 0 for j ~ 3, applying Lemma 2 we get that M(pCI» is the parallelotope P(t(1), given by

bP»

1'1

-

. 2 ,... ,n }=

1'1

[

.

1'1 Hence, by Theorem I, it turns out that

tCO) _ 1'j tCO) ) 1'1 I

;(0) = !!L,

r Si So Sj

J.l(So n S2) = -J.l(SI n S2) .

~1) = O~O) + .!-t~O)(co _ poO~O» .

t<1)

=

Sj.

From (15) we

- 1'1 1

t )~l)

if i :f i· ,j* if 1'i. > r otherwise if 1'j. > 1 otherwise

Si

J.l(SonSd = -J.l(SI nS2)

t(l) - .!-t(O) 1

Let

=

Then, we have the following result.

P

[co] C2

thus proving the Lemma.

PI M and where PI and P2 are such that PI P2 M , respectively. Finally, define Cl Cl, C2 C2 and Co co·

P2

0(0) _ ...,,-l ti O)(co _ poO(O» = c /1 1 c

1'j .

Clearly, by construction So is oriented and So = So whenever So is oriented. Moreover, since Po E span {ml' m2} there exists a vector Po such that Po = PoM. From (19), it turns out that

=

-

M(V)

j=3

Po

c

=

n

r = 1+ L

#.1) _

~O) _

1 =

] 'Y2 t 0 )

1'1

~O

[ _ ] -1

P2

for each oriented parallelotope P containing V. If i* 0, then P* M(pU» for each oriented parallelotope pU) = {ni=O,i;tjSd. Thus, the result follows by applying Theorem 1 in the parameter space, with the additional constraints that

=

519

=

M(P*) = P*. If i* :f. 0, then define the oriented parallelotopes

Gu, G., and P. Khargonekar (1992). Linear and nonlinear algorithm; for identification in H 00 with error bounds. IEEE Trans. on Automatic Control, 37,953-963.

= So n {n?=O,i~io,j;fjSd . Once again, M(PU» = P* and the result follows PU)

from Theorem 1.

Helmicki, A. J., C. A. Jacobson and C. N. Nett (1991), Control oriented system identification: a worst case/deterministic approach in H oo . IEEE Trans. Automat. Contr., 36, 1163-1176.

•

4. COMMENTS AND CONCLUSIONS Lemma 1 and Theorem 2 provide basic tools for constructing an algorithm which recursively updates the outer parallelotopic bound for the feasible parameter set. The result of Theorem 2 guarantees that at each iteration, the updated estimate is the parallelotope whose value set at a given frequency is of minimum area, within the class of oriented parallelotopes. The obtained parallelotope is of minimal volume within the same class. Also, it is conjectured that the parallelotopic estimate minimizes the value set area even if the class of admissible parallelotopes is the set of all parallelotopes.

Kosut, R. L., M. K. Lau and S. P. Boyd (1992). Set-membership identification of systems with parametric and non parametric uncertainty. IEEE 1rans. on A utomat. Control, 37, 929-941. Milanese, M., and A. Vicino (1993). Information Based Complexity and non parametric worstcase system identification. Journal of Complexity, 9, 427-446. Vicino, A., and G. Zappa (1993). Sequential approximation of parameter sets for identification with parametric and non parametric uncertainty. Proc. 9£-nd Conference on Decision and Control, pp. 2044-2049, San Antonio, USA.

It is clear from Theorem 2 that the obtained so-

lution represents a trade-off between the value set area and parallelotope volume reduction. More specifically, if i* = 0, i.e., if both 1'1, 1'2 < r, then the algorithm does not provide any reduction in the area of the model value set, while the volume reduction rate is provided by 1/1'jo. However, if either 1'1 or 1'2 > 1'jO , a greater volume reduction could be obtained at the expense of an area increase in the frequency domain by replacing S1 (or S2) with So. If i* :f. 0, the area reduction rate is provided by r hio and the volume reduction by r / (-Yi 1'j 0). Since r > 1'j 0, a greater volume reduction rate 1/1';0 with the same area reduction could be obtained by the non oriented parallelotope p(i") = {ni=O,j;fioSd.

Wahlberg, 8., and L. Ljung (1992). Hard frequency-domain model error bounds from least-squares like identification techniques. IEEE 1h1ns. on A utomatic Control, 37, 900912. Younce, R. C., and C. E. Rohrs (1990). Identification with non-parametric uncertainty. Proc. £9th CDC, pp. 3154-3161, Honolulu, USA.

o

As a last observation, notice that, if u(k) = cos (wk +4» and a FIR parameterization is adopted, the regressor 4>'(k) E span {m1(w), m2(w)}. Hence, the strip S(k) is oriented and, at any step, optimal area reduction implies optimal volume reduction. For a general input signal, it is believed that using the algorithm to minimize the area of the value set at different frequencies could lead to a consistent frequency domain uncertainty reduction.

REFERENCES Chen, J., C.N. Nett and M.K. Fan (1992). Optimal non parametric system identification from arbitrary corrupt time series: a control oriented approach. Proc. 199£ Amer. Contr. Conf., pp. 279-286. Dahleh, M., and R. Smith, Eds. (1993). The Modeling of Uncertainty in Control Systems. Springer Verlag, London. 520