Robust Parameter Estimation for Linear Models with Set-membership Uncertainty

Copyright © IF AC Identification and System Parameter Estimation , Beijing, PRC 1988

ROBUST PARAMETER ESTIMATION FOR LINEAR MODELS WITH SET-MEMBERSHIP UNCERTAINTY G. Belforte*. R. Tempo** and A. Vicino* *Dipartimento di Automatica e Informatica, Politeenico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy **CENS-CNR, Politeenico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy

Abstract Y A A

In this paper we address th .. problem of parameter estimation of a linear model P where the input matrix A is known and the additive uncertainty P is assumed to be unknown but bounded in an I"" norm by a given constant (. In this case for given data y the set of all admissible parameters A consistent with the given model equation and the bound { is a convex polytope, In this paper we attempt to estimate optimally, in the sense of maximizing the volume, a rectangular box which is guaranteed to lie completely inside the true region, The inner bounding region obtained in this way can be used together with an outer bounding region, in order to approximate as closely as possible the boundary of the actual admissible parameter polytope, A computable algorithm, using linear programming, is proposed and some numerical results are given to illustrate the behavior of the technique.

+

Key Words Bounding Set .

Robust Parameter Estimation, Linear Models, Set-Membership Uncertainty, Inner

1. Introduction (3) Many different problems of theoretical importance and of practical interest can be formulated in terms of a linear relation between an unknown parameter vector A = [AI, ... , An)T in an n-dimensional space, and a vector Y = [YI,"" Ymf of known quantities in a m -dimensional space Y=

AA.

where { represents a known quantity. A case of great concern is when I"" norms are adopted; in this case an Unknown But Bounded Error (UBBE) assumption is used in the sense that each component of the error vector is known to be bounded by given values. Motivation for this kind of error representation is the fact that, in many practica.l cases, the UBBE information is closer to the actually available information on the measurement error. As an example consider the case in which measurements are obtained by means of instrumentation whose precision is defined in class, or the case of truncation (A I D converters) and roundoff errors, or

(1 )

In the following we will call Y the measurement vector and we will mainly refer to parameter identification problems, for which A represents the system model. However we recall that many other problems such as state estimation and time series forecasting in systems and control area (Milanese and Tempo, 1985; Milanese and Belforte, 1982), indirect measurements in measurements area (Smit , 1983) , can be described by relation (1). In this paper we will consider the case in which m 2: n; in such conditions an exact solution of (1) cannot be found in most practical applications owing to the perturbed available information that is caused by, or is assumed to be adequately represented by an additive measurement error on y. Then relation (1) must be rewritten as: Y=

AA + P

the case of tolerances that are of paramount importance in mechanical applications. In this context the identification process consists in finding the set of all the parameter values consistent with the system model and with the error model (parameter uncertainty set). Any element of the set represents a possible estimate although the center or the minimum norm element of the set deserves interesting optimality properties (Traub and Wozniakowski , 1980; Micchelli and Rivlin,1977; Kacewicz et aI. , 1986). The size of the set represents a measure of the estimates reliability. Unfortunately, an exact representation of the parameter uncertainty set is in general not simple. In fact, this set is a convex polytope in the parameter space whose boundary is defined by a set of 2m inequalities deriving from subsitution of (2) in (3). It is therefore convenient to look for simpler although approximate descriptions of the parameter uncertainty set. To this extent the use of simple shaped sets like ellypsoids or orthotopes bounding the parameter uncertainty set has been proposed (Milanese and Belforte, 1982; Fogel and Huang , 1982). The computation of the orthotopic bound to the parameter uncertainty set is exactly achieved solving 2n linear programming problems to com-

(2)

where P = [PI,' .. ,Pm]T represents the m-dimensional error vector . As far as some information about the error model is not available no solution to relation (2) can be found. The classical approach to the problem is to assume that the error vector P is statistically modelled by a known (at least partially) distribution. In more recent years a new approach, referred to as set membership error description, has been investigated. In this case the error vector P is assumed to be the element of a given set; a special case of interest is when the shape of this admissible error set is described in an Lp norm

809

G. Belforte, R. Tempo and A. Vicino

810

pute the maximal and minimal values of each parameter. In this way the parameter uncertainty intervals obtained can be regarded as parameters confidence intervals with a confidence level equal to one. This kind of characterization of the parameter uncertainty set is of great importance in many identification applications when an estimate and its variability are of interest; however, the information derived from the bounding orthotopes is obviously incomplete and other pieces of information could be sought. In particular, inner orthotopes, i.e., boxes completely contained in the parameter uncertainty set, could give important additive information of primary importance for some applications. As an example, consider the case in which the identification goal is the determination of the parameter values consistent with given specifications; in this case an inner box, whose cent er and edges can be regarded as the nominal parameter values and parameters admissible tolerances, is more useful than a bounding orthotope whose elements can lie outside the parameter uncertainty set. Indeed, all the elements of the inner box belong to the parameter uncertainty set and can be regarded as a robu&t e&timate. Aim of this paper is to discuss the problem of determining inner orthotopes of maximal volume. Some algorithms for their computation are introduced and, in particular, it is shown that the computational complexity is related to some possible additive information about the location of the centers of the orthotopes and / or their shape.

Ily - A).II::;'

~

(4)

f

where f is some fixed positive number. An algorithm

Z i.e., it provides an approximation
Ily - A).I I::;'

T(y) = {>. EA:

~

f}

(5)

We assume that T(y) is nonempty, i.e. , that the model structure is able to represent all the data y belonging to a set Yo <; Y. The set T(y) is a convex polytope and represents the set of all admissible parameters). compatible with the information A, the data y and the bound on the noise

E-

Now we present a simple parameter estimation problem showing the role of the various spaces and operators introduced above. Example 1 Assume that the space A is a two-dimensional space of real time-functions

A

= {f(t)

: f(t)

= ).1 + ).2t,

t

2': O}

(6)

This section provides formal definitions and notations used in the paper. Let A be a linear normed finite dimensional space over the real field. Consider a given linear operator S, called a wlution operator, which maps A into Z

Notice that A can be equivalently identified either with the space of real time functions f(t) or with the space of the parameters). defining f(t) . We consider the problem of estimating the parameters ). knowing three perturbed measurements of f(t) at times t = {O, 1, 2} with some norm bounded uncertainty. In this case the solution space Z coincides with the space of the parameters A to be estimated and Y is a three-dimensional measurement space

S:A -> Z

A = Z == R2 and Y = R3

2. Definitions and Notation

where Z is a linear normed finite dimensional space over the real field . Our aim is to estimate an element S). of the space Z knowing approximate information about the element ). . Define a linear operator A, called information operator, which maps A into a linear normed finite dimensional space

Thus the solution operator S : R2 operator

A : A -> Y. In this paper we assume that dimY 2': dimA and that the information operator A is a complete information , i.e., it is a one-to-one mapping. The former assumption is satisfied in many practical problems we are interested in , i.e., system parameter estimation and time series prediction. The latter assumption means that we can exactly recover). from the knowledge of the exact information A),. In the following we assume that A and Z are equipped with 1= norms and the space Y is equipped with an I: norml. The use of an norm enables us to consider different bounds on every measurement . In general, in the presence of noise , exact information A), about). is not available and only perturbed information y is given. In this context we assume that the uncertainty of information p = y - A). is unknown but bounded

I:

1 The

l~

norm is defined as

R2

IS

the identity

(7) and the information operator A is the sampling operator defined by

Af( t) = [J(O),J( 1), f(2}f =

Y

->

[).I,).I

+ ).2 , ).' + 2).2]T (8)

The information Af(t) is corrupted by an additive noise p == y - Af(t) bounded by

with weights given by

WI

==

W2

Ilpll:

~

==

= 1. The available data y are

W3

y ==

f

= 3

[4, 1, 6]T

(9)

(10)

It results that the outer box bounding the region T(y) is defined by the following inequalities

(11) while the maximal inner box with squared shape is gi ven by

(12) A geometrical description of inner and outer boxes is presented in Figure l. In this case it can be observed that the description of outer bounding box is not tight to the true region and the computation of inner bounding box may

Robust Parameter Estimation for Linear Models be suitable. Several different algorithms may be chosen to provide specific estimates of the parameters >. as, for example, the minimum norm >.MN '" 13.5,0.1 or the central estimate >.C '" [3.,1.1. The general formulation previously introduced can be used for problems different than estimation of the parameters of a dynamic system. For example, the use of this approach may be useful for time series prediction (Vicino et aI., 1987). Knowing approximately m past samples of a time function f( t)

[y(tl ),y(t 2 ), ... , y(tmW '" .4>'

Notice that in (17) the values of the weights Vi determine the shape of the box. The following two problems are of interest in robust estimation contexts . • p.roblem 1 Given an admiJJible vector >.0 (i_e., >.0 E T(y) ), find the maximal box of given shape centered at >.0 and entirely contained in T(y). • Problem 2 Find the maximal box of given shape entirely contained in T(y).

+ p '"

'" [f(td.J(t 2 ), ... .J(tmW + P

811

(13)

our goal is to forecast f(t) at times {t m+1 ,t m+2 , •.• ,tm+r}. It results that the solution operator 5 is the shift operator defined by

In this case an inner box estimate provides upper and lower bounds of the function f(t) over the whole forecasting horizon {tm+ 1, . . . , tm+r}'

Now we show that Problem 1 admits a one shot solution and Problem 2 admits a solution in the form of an LP problem. Problem 1 Let us translate the origin of the space A in >.0 and define the coordinates with respect to the new origin as

(18) The polytope T(y) in the new coordinate axes becomes

3. Robust Parameter Estimation of Linear Systems As previously mentioned, when (2) and (4) hold the region of admissible parameters T(y) given in (5) for any given measurement vector y E Yo is a convex polytope. It has been shown in (Milanese and Belforte, 1982) that the uncertainty intervals for each parameter >'i compatible with the available information can be computed by solving suitable Linear Programming (LP) problems. More precisely, if we denote by >'i' Pt1) the minimum (maximum) value of >'i in the set T(y) (i '" 1, ... ,n), we know that these extremal values are given by

>';" '" inf >'i

(15)

>.t' '" sup '\;

(16)

To(y) '"

P. E A: Il y -

A~ - A>.°l l:' S

"' {~ E A:y;-S ta;j:i.j : ; yt,

i "' I, ...

}= l

f} '"

,m}

(19)

where aij are the entries of the matrix A and n

Yi - Laij).~

+ Wif, i = 1, . .. ,m

(20)

i = 1, .. . ,m

(21)

;=1 n

Yi - 2:aijAj -

Wi€,

; =1

The solution to Problem 1 for a given measurement vector y can be obtained by solving the optimization problem

i. ET(y)

sup I I~II:;"

(15) and (16) are LP problems with n parameters and 2m inequality constraints. It can be observed that the computation of >'i' and >.t! leads to the estimation of the minimal box bounding the polytope T(y) from outside. Moreover, the geometrical cent er of this outer box is the center ( in the Chebychev sense) of the set T(y) when an I"" norm in the parameter space is used and can be considered as a point estimate of the parameter vector >.. This estimate enjoys strong optimality properties as shown in (Milanese and Tempo, 1985; Kacewicz et aI., 1986). Nevertheless, the outer box estimate, though providing interesting information on T(y) , cannot be considered a robuJt set estimate since, in general , it may contain also points not belonging to T(y). On the contrary, a box completely contai ned in T(y) can be considered a robuJt estimate of the parameter uncertainty set. In the following we attempt to compute this robust estimate imposing at the same time the maximality of the inner box. We assume that the box looked for belongs to a family of boxes B~(6) of fixed shape and centered at some point >" B~(6)"'{>. E A:ii>. - >"I I :;" S 6},

6 ::: 0,

V l , " " V n ::: O

(17)

(22)

.\ET,(y)

>' ET(y)

which is equivalent to sup 6> 0

(23)

II ~ II ;;" S 6 ~ E To(y) Observe that, due to the convexity of To(y) , one for each oC the double sided inequality constraints in (23) are necessarily sharp in correspondence to the solution of the infimum problem (23). The tangency point of the box II~II;;" ::; 6 with the boundary 8To(y) lies on one of the hyperplanes defin ed in (19). Consider th e two hyperplanes determined by the i-th measurement n

2: a ;)j

y;+

(24)

y;-

(25)

i=1 n

2: a ;)j j=1

By substituting the active constraint equat ions of (23) in (24) and (25), we obtain

C. Belforte. R. Tempo and A. Vicino

812

Setting (26)

v

n

o

L

(27)

la;;v;1

IVI,'" IYI

y

[Yl -

;=1

la;;v;1

(28)

la;;v;1

(29)

n

Iy;-I I L

fWl, . ..

(39) ,ym

+ fWm]T

,Ym -

(40) (41)

fWm]T

and using (36) and (37) problem (33) reduces to

n

lytl / L

,vnf

+ fWI,'"

y+

;=1

and hence

(38)

{la;;I}

IAI

i =l

;=1

y-

+

(42)

sup 0> 0 y - S AA c S y+ olAlv S AA c S y+ - olAlv

where ot and 0;- are values of 0 obtained by solving (26) and (27) respectively. Now we state a theorem following from the above considerations. Theorem 1 Let ot,o;- be given by (28) and (29), for i = 1, ... ,m. The maximal inner box centered at AO iJ given by

An easy inspection of (42) shows that the second set of constraints is redundant , so that the original problem can be reduced to

(30)

which is an LP problem with (n + 1) parameters and (2m + 1) inequality constraints. The above results can be summarized in the following Theorem 2.

where 0" = min {

min

; E{I •...• m}

ot,

min

; E{I , ... ,m}

o;- }

(31)

sup 6> 0 y - + 61Alv S AA c S y+ - 61Alv

(43)

Theorem 2 Let li and Xc Jolve the LP problem (4:1) with (n + 1) parameterJ and (2m + 1) conJtraintJ , Then, the maximal box B~(6) entirely contained in the uncertainty pa-

Problem 2 To find the maximal box of given shape entirely contained in T(y) with variable center AC , we will consider again the class of boxes (17). Given a measurement vector y, the solution can be obtained by solving the following maximization problem sup AC E T(y) A E T(y)

(32)

which is equivalent to (33)

sup 0> 0

IIA - AC II;;., S o lIy - AAII: S t lIy - AAclI: S ( (33) is an LP problem in (2n+1) parameters and (4m+2n + 1) inequality constraints. Now we show that the complexity of this problem can be reduced significantly. In particular , we observe that, as in Problem 1, the values of the parameters for which the infimum (33) is achieved are such that n out of the 2n inequality constraints in the second set of constraints in (33) are sharp. Consider the generic i-th constraint of the third set in (33) and the constraints included in the second set

rameter Jet T(y) iJ 9iven by

( 44)

Remark It is worth while to observe that the solution presented above allows to solve , as a special case, an additional problem which is generally of concern in gaining information about the parameter uncertainty set. This problem consists in evaluating some diameterJ of the polytope T(y), called coordinate diameterJ , consisting in the maximum distances among points on the boundary 8T(y) along the coordinate directions of the space A. We can compute the coordinate diameters by solving again problem (33) , setting all weights v; to zero, except for that relative to the parameter of the direction considered , Thus, in order to compute the diameter along the direction A;,j E {I , . . . , n} we have to solve the problem sup 6>/ . {

, .

y-

6> 0 + oIA; lv; S A.A c :s

6

(45)

y + - oIA; lv;

where IA;I represents the j -th column of matrix IAI . If 0; and AC are the values of the parameters at the supremum in (45), the coordinate diameter looked for is given by

n

y; -

tW;

S

L a;; oX; S y; + tw;

(34)

j=1

oX;

= Aj ± ov;, j = l , ... , n

(35)

By substituting (35) in (34) we obtain two independent constraints for each i E {1, . .. , m}

a;; IAj

+ Dv; sgn( a;;) ] S

y;

+ (lL';

(36)

;=1 n

L ; =1

j}

(46)

In the next section we present a numerical example of application of the above results to a dynamical system.

4. Numerical Example

n

L

Ai = X~ , i = 1, ... ,n,i ;i

a;;IAj - Dv; sgn(a;;) ] :::- Yi -

fll'i

(37)

In this section we present the computation of inner and outer boxes and the coordinate diameters of T(y) for the following MA(3) model studied in (Milanese and Belforte, 1982)

813

Robust Parameter Estimation for Linear Models

y(k) = AIU(k) + A2U(k - 1) + A3U(k - 2)+e(k) k = 1, . . . ,50 (47) where the input u(k) is a known sequence generat ed by a gaussian distribution and e(k) is a sequence of uniformly distributed measurement noise. The following nominal values of parameters are considered Al = 1.25

A2

=

References [1] E. Fogel and F. Huang, "On the value of information in system identification-Bounded noise case," Automat. ica, vol. 18 , pp. 140-142 , 1982. [2] B.Z. Kacewicz , M. Milanese, R . Tempo and A. Vicino, "Optimality of Central and Projection Algorithm," Sy~tem~ and Control Letter~, vo!. 8, pp. 161·171 , 1986.

2.35

A3 = 0.50

[3] C .A. Micchelli and T .J . Rivlin, " A survey of opti · mal recovery", in : C.A. MiccheUi and T .J . Rivlin , Eds. , Optimal E~timation in Approximation Theory (Plenum, New York, 1977) pp. 1-54.

Table 1 presents the computation of the unceratinty intervals of the estimated parameters of the model (47) for f

= .5 ~l

I

~3

~2

i~ner box 1 [1.254 ,1.263]1 [2.341 ,2.350] outer box [1.202,1.288] [2.332 ,2.367]

[4] M. Milanese and G . Belforte, "Estimation theory and uncertainty interval evaluation in presence of unknown but bounded errors: Linear families of models and estimators ," IEEE Tran~. on Automatic Control, vo!. AC27, pp. 140·142, 1982.

[0.494 ,0.503] ! [0 .461 ,0.513]

Table 1

[5] M. Milanese and R . Tempo, "Optimal Algorithms Theory for Robust Estimation and Prediction" , IEEE Tran~ . on Automatic Control, vo!. AC·30 , pp. 730-738, 1985.

Table 2 shows the coordinate diameters of the uncer tainty parameter set. ~

Dl

D3

[ diameters : [1.254,1.281 ] [ i2.333,2.36Z j [ [0.494 :0.512]

I

[6] M . K . Smit , " A novel approach to the solution of indi rect measurement problems with minimal error propagation " , Mea~urement , vo!. 1, pp. 181·190 , 1983.

Table 2

Acknowledgements

[7] J .F. Traub and H. Woiniakowski , A General Theory of Optimal Algorithm ~ (Academic , New York, 1980).

This work has been partially supported by funds of Ministero della Pubblica Istru zione and of CENS-CNR. The third author (A . Vicino) is presently with the Dipartimento di Si~temi e Informatica, Univer~itd di Firenze, Firen ze (Italy) .

[8] A. Vicino , R. Tempo, R. Genesio and M. Milanese, "Optimal Error and GMDH Predictors: a Comparison with some Statistical Techniques ," Th e International Journal of Foreca~ting, pp. 313-328 , 1987.

5

3

--------1 I

2

)...C -5 -4 -3

-2 -1

I I

"MN: A1 5

FIGURE 1 . Inn.r end out.r box •• for Ellempl. 1.