Linear Time-Invariant Models of Non-Linear Time-Varying Systems

Time-varying Systems

Lennart Ljung

Department of ElectricalEngineering

Linkoping University, SE-581 83Linkoping, Sweden

WWW: http://www.control.isy.l iu. se

Email: ljung@isy.liu.se

October2, 2001

REG

LERTEKNIK

AUTO

_{MATIC CONTR}

OL

LINKÖPING

Report no.: LiTH-ISY-R-2363

Forthe 2001 European ControlConference, Porto, Portugal. Also

published in The European Journalof Control, Vol7, issue 2-3,

2001

Time-Varying Systems

LennartLjung

Div. of Automatic Control

LinkopingUniversity

SE-58183 Linkoping,Sweden

email: ljung@isy.liu.se

April 10, 2001

Abstract

Thestandardmachinery forsystemidentication of

lineartimeinvariant(LTI)modelsdeliversanominal

model and acondence (uncertainty)regionaround

it,basedon(secondordermoment)residualanalysis

and covarianceestimation. In most casesthis gives

anuncertaintyregionthattendstozeroasmoreand

moredatabecomeavailable,evenifthetruesystem

isnon-linearand/ortime-varying. Inthispaper,the

reasonsforthisaredisplayed,andacharacterization

of the limit LTI model is given under quitegeneral

conditions. Various ways are discussed, and tested,

to obtain a more realistic limiting model, with

un-certainty. These should re ect the distance to the

true possibly non-linear, time-varying system, and

also form areliablebasis forrobust LTI control

de-sign.

1 Introduction: The Fiction of

an LTI System

Linear,Time-invariant(LTI)descriptionsof

dynam-icalsystemsareclearlythebreadandbutterof

con-troltheory. Nevertheless,theyarestill action: No

real-life systemis exactlylinear and time-invariant.

So,althoughtherearenoLTIsystemsoutthere,LTI

modelsasabasisforcontroldesignhaveprovedtobe

of enormousvalue. There are basically tworeasons

for this: (1) anLTI model may be agood

approxi-mationofareallife systemand(2) feedbackcontrol

is forgiving, in the sense that youcan achievegood

controlbasedonquiteanapproximatemodel.

SystemIdentication oersaneÆcient machinery

to estimate LTI models from observedinput-output

data. Thismachinerywillbebrie ysurveyedin

Sec-tion 2. Identication techniques deliver a nominal

LTImodel,withanassociateduncertaintyregion,

re- ectingtheestimatedstatisticalcondenceregionof

theestimatedparameters. Itis intuitiveto visualize

the delivered model as a band around the Nyquist

curve or as bands in the Bode plots. These

con-denceregions are deemed to bereliable (or at least

"not falsied") if certain model validation tests are

passed. A typicalsuch test is to check the

correla-tionbetweenthe modelresiduals (prediction errors)

andpastinputs,aswellasthecorrelationamongthe

modelresidualsthemselves.

It is important to realize that the LTI

identica-tionmachineryisalwaysabletodeliveranunfalsied

linear model with decreasing uncertainty regionsas

moreand moredatabecome available,regardlessof

the character of the system. The reason for this is

thatLTI-techniques(i.e. secondordermoment

tech-niques)cannotdistinguishasystemfromitsLTI

sec-ondorderequivalent. Thedetailsofthisarediscussed

Example 1.1: Rotation ofa Rigid Body

Consider the rotation of a rigid body around a xed

point. Theinputu isamomentappliedalongacertain

axis. Theoutputyistheangularvelocityaroundoneof

theprincipalaxesofthemomentofinertia. Ifuisapplied

aroundthesameprincipalaxis,therotationisdecoupled

andthereisalinearrelationshipbetweenuandy.

How-ever, if u is appliedalong another axis, thesystem will

notbelinear,unlessthebodyissphericallysymmetrical.

The reason is that there are non-linear cross-couplings

betweentheprincipalaxesofthemomentofinertia. To

be specic, weconsider athinand long bodywith

mo-ments ofinertia being 0.11, 100.01, and 100.10,

respec-tively, along the principal axes. There is also aviscous

dampingfactorof0.01aroundeachoftheaxes. We

per-formtwoexperiments:

 A.Theinputmomentisappliedaroundanaxisthat

deviatesfromtheoutputprincipalaxisby0.003

ra-dians.

 B.Theinputmomentisappliedalongtheline that

deviatesfromtheoutputprincipalaxis by0.1

radi-ans.

Therstsystemwouldthenbe"ratherlinear",whilethe

secondoneishighlynon-linear.

Ineachcasealowfrequencyrandominputwaschosen

and10000datapointscollected. Portionsofthedataare

showninFigures1-2. Thelinearidenticationprocess

selected a thirdorder BJmodel inboth cases. Results

fromresidualanalysisinthetwocasesareshownin

Fig-ures 3 -4. Neither give any reason to reject the

mod-els. AmplitudeBodeplotswithcondenceregions

corre-sponding to3 standarddeviations are shown inFigures

5 -6. Thedelivered picture is clear: Bothsystemscan

condentlybedescribedbyLTImodelswithonlyminor

uncertainty.

ForexperimentA,themodelisabletodescribe100%

oftheoutputvariationbyone-stepaheadpredictionand

99.39%inapuresimulation. Thecorrespondinggures

forexperimentBis99.97%and0.13%,respectively.The

last gure is low, but theLTI-identication machinery

explainsthedeviationasnoisethatis uncorrelatedwith

the input and can be described as ltered white noise

disturbances, giving rise to some limited uncertainty in

theestimatedmodel.

Ifweknowthatthedatacollectionhasbeenessentially

noise-free,thisinterpretationshouldcausesomeconcern,

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

3.8

4

−0.02

−0.015

−0.01

−0.005

0

y1

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

3.8

4

−5

0

5

10

u1

Figure 1: Portions of thedata. ExperimentA. The

upperplotisoutput(angularvelocityaroundthe

sec-ond principal axis.) andthe lowerplot isthe input

(appliedmomentaroundanotheraxis.)

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

3.8

4

0.006

0.008

0.01

0.012

0.014

0.016

0.018

0.02

y1

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

3.8

4

−10

−5

0

5

10

u1

Figure2: Portionsof thedata. ExperimentB.

0

5

10

15

20

25

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Correlation function of residuals. Output y1

lag

−25

−20

−15

−10

−5

0

5

10

15

20

25

−0.03

−0.02

−0.01

0

0.01

0.02

0.03

Cross corr. function between input u1 and residuals from output y1

lag

Figure 3: Result of residual analysis. Experiment

A.Theuppercurveshowstheautocorrelationofthe

residuals. Thelowerplotshowsthecrosscorrelation

betweenresidualsandinputs. Theshadedzoneisthe

condenceregion.

0

5

10

15

20

25

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Correlation function of residuals. Output y1

lag

−25

−20

−15

−10

−5

0

5

10

15

20

25

−0.03

−0.02

−0.01

0

0.01

0.02

0.03

Cross corr. function between input u1 and residuals from output y1

lag

Figure4: Resultofresidualanalysis. ExperimentB.

Sameexplanationasin Figure3.

10

0

10

1

10

2

10

−5

10

−4

10

−3

10

−2

Frequency (rad/s)

Amplitude

From u1 to y1

Figure 5: Amplitude Bode plot of the obtained

model, with condence regions corresponding to 3

standarddeviationsmarked. ExperimentA.

10

0

10

1

10

2

10

−5

10

−4

10

−3

10

−2

Frequency (rad/s)

Amplitude

From u1 to y1

Figure 6: Amplitude Bode plot of the obtained

model, with condence regions corresponding to 3

did the simulation and knowthat the system incase B

is highly non-linear, the Bode plotin Figure 6 and the

residualtestplotinFigure4shouldcallforfundamental

concern.

Theexamplepointsoutafundamental

shortcom-ingof thestandardLTIidentication process: With

increasing amounts of data, models will be

deliv-ered with uncertainty zones converging to zero in

Nyquist/Bode diagrams. This does not rhyme well

with our knowledge that while LTI models may be

good approximations, no real life system is exactly

LTI.Itwouldbemuchmoresatisfactoryifthe

deliv-eredLTI model has someremaininguncertainty, no

matterhowmanydataitisestimatedfrom.

Thetopicofthis contributionistodiscuss this

is-sue.

Dealingwithremainingbiaserrorsinmodelsisby

nomeansanewproblem. There aremany

contribu-tionsin theliteraturethatdealwiththeproblem to

livewithbothbiaserrorsandtheclassicalstatistical

varianceerrors. Wecouldpointto, amongmany

ref-erences,[1]foracharacterizationofthebiaserrorin

thefrequency domain, [2]and [3] fortheconcept of

stochastic embedding, [4] for model error models, [5]

fortotalerrorestimates,[6],[7]and[8]formore

de-terministicmeasures, [9]forexplicit analysis ofbias

andvariancecontributions,[10]formodel

approxima-tionstailored to control design,and [11] forexplicit

robustnessmeasures foridentied models.

Most of these references, however, deal with the

problem that the model is of lower order than the

truesystem,whichstillisassumedtobegivenasan

LTIdescription. Inthispaperwewillspecically

dis-cussmodeldiscrepanciesthat arecaused bysystems

that are more diÆcult to describe. An early

treat-ment ofLTI modelsand ill-dened systemsis given

inChapter 8of[12].

LTI Models

A generalLTI-model of adynamical systemcan

al-waysbedescribedas

y(t)=G(q;)u(t)+H(q;)e(t) (1)

Here, q is the shift operator, and G and H are the

transfermatricesfromthemeasuredinputuandthe

noisesource e, which is modeled aswhitenoise

(se-quence ofindependentrandomvariables). For

nota-tionalconveniencewewillfromnowononlyconsider

Single-Input-Single-Output systems, but the theory

isthesamein themultivariablecase.

Weshallalsousethefollowingshorthandnotation

forthecorrespondingfrequencyfunction

G

 =G(e

i!

;) (2)

Thetransferfunctionsareparameterizedbya

nite-dimensionalparametervector,andthis

parameter-izationcanbequitearbitrary. Forblack-boxmodels,

itiscommontoparameterizeGandHintermsofthe

coeÆcientsofnumeratoranddenominator

polynomi-als,perhapsconstrainingGandH tohavethesame

denominators. This leadsto well establishedmodel

classes,knownundernameslikeARX,ARMAX,OE,

BJ,etc.

Anotherpossibilityistoparameterizethemodelas

astate-spacemodel, indiscreteorcontinuoustime:

x(t+1)=A()x(t)+B()u(t)+K()e(t) (3a)

y(t)=C()x(t)+D()u(t)+e(t) (3b)

which gives G(q;)=C()(qI A()) 1 B()+D() H(q;)=C()(qI A()) 1 K()+I

Theparameterizationofthestate-spacematricescan

beofblack-boxcharacterascanonicalforms,oreven

llingallthematriceswithparameters. Itcanalsobe

in termsofagrey-box,where physicalinsight

(typi-callyin continuoustimedescriptions)isused,mixed

estimate the parameters in (1) based on observed

input-output sequencesfy(t);u(t);t =1;2;:::;Ng.

Amongmanysuggestedalgorithmsforthis,two

ma-jorapproachesaredominatingtoday:

 Sub-spacemethods

 Predictionerrormethods

Sub-space methods, e.g. [13], [14], [15], canbe

de-scribedasrstestimatingthestatesequencein(3a)

and then treatingthe two equations, with assumed

knownx,as linearregressionstondthestatespace

matrices.

Prediction errormethods rst determinethe

pre-dictionerrorsassociatedwith (1):

"(t;)=H 1

(q;)(y(t) G(q;)u(t)) (4)

ThisrequiresbeconnedtoaregionD,sothatthe

lters H 1

and H 1

G are stable. Then the  that

minimizesthenormoftheerrors

^  N =argmin 2D V N () (5a) V N ()= 1 N N X t=1 " 2 (t;) (5b)

isdetermined,typicallybynumericalsearch. Agood

combinationofthetwoapproaches,in theblack-box

case, is to initialize the search at the estimate

pro-videdbythesubspacemethod.

How will these methods perform? Well, that

de-pendsontheinput-outputdata. A typicalapproach

to analysis is to assume that the data indeed have

been generated by a system like (1) for some

par-ticular parameter vector 

0

, and for e being a

se-quence of independent random variables. In that

case the asymptotic statistical properties

(conver-genceandasymptoticdistribution)of ^



N

canbe

cal-culated readily. We refer to [16] for a

comprehen-sive analysis of this kind, as well as for more

de-tailsonmodelstructuresandestimation techniques.

Justonethingwillbepointedout,though: Itispart

ofthestandardLTI-identicationmachineryto

com-putetheresultingresiduals:

"(t)="(t; ^



N

) (6)

pastinputsu(s);st andif theyare mutually

un-correlated. If such a residual analysis testis passed

(i.e. there is no convincingstatistical evidence that

correlationispresent),theassumptionofatrue

sys-tem within (1) corresponding to a particular value



0

is \notfalsied", and thedistribution of ^



N 

0

can be calculated using the aforementioned theory.

Thismeansthatacondenceregionforthetrue

sys-temcanbeestimated. ThedeliveredLTImodelthus

comeswithaqualitytag,correspondingtocondence

regionsaroundtheestimate. This wasdepicted,e.g.

inFigure5.

Instead of reviewing this standard material, we

shall in this paper develop an independent analysis

of thelimitof theprediction errormethod estimate

^



N

. This willuseminimalassumptionsonthe

prop-ertiesoftheinput-outputdata. Inparticular, itwill

not be assumed that they have been generated by

an LTI system, and it will not employ a stochastic

framework. Some related results were presented in

[17].

3 Second Order Equivalent LTI

Models

3.1 Quasistationary Signals

Adeterministicsignalz(t)willbecalled

quasistation-ary,[16],if

jz(t)jC ;8t forsomeC<1 (7a)

lim N!1 1 N N X t=1 z(t)z T (t )=R z (); exists8 (7b) IfR z

issuchthattheZ-transform

 z (z)= 1 X = 1 R z ()z  (8)

is well dened on the unit circle, we call 

z (e

i!

)

the spectrumor spectral density of z. 

z

(z) will be

called the spectral function. It can be shown that

R

z and

z

possessallthe properties normally

forstationarystochasticprocesses. InSection4.2we

shallspecicallyprovehowtheytransformunder

lin-earltering.

We will also use the following standardconcepts:

Alter G(z)= 1 X k = 1 g k z k (9) willbecalled  stableif P jg k j<1  causalifg k =0; k<0  strictlycausalifg k =0; k0  anti-causalifg k =0; k>0.

Moreover,afamilyoflters

G  (z)= 1 X k = 1 g  k z k ; 2D (10)

iscalleduniformlystable if

1 X k = 1 sup 2D jg  k j<1 (11)

3.2 Description of Systems that

Pro-duce Quasistationary Data

Lettheinput-outputdatacollectedfromtheprocess

befu(t);y(t);t=1;2;:::g. Let

z(t)= 

y(t)

u(t) 

Assume that the data are quasistationary and that

thespectralfunction

 z (z)=   y (z)  yu (z)  uy (z)  u (z)  (12) iswelldened.

Now,dospectralfactorization



z

(z)=L(z)L T

(1=z)

sothatL(z)andL (z)arestableandcausal2-by-2

transferfunction matrices. Then dene

P(z)=   yu (z)  y (z)  L T (1=z) 1 = 0 X k = 1 p k z k + 1 X k =1 p k z k =P (z)+P + (z) whereP +

(z)isthestrictlycausalpartofthelefthand

side. NextdeneW

u andW y by P + (z)L 1 (z)=  W u (z) W y (z)  (13) Byconstruction W u and W y

will be strictly causal,

i.e. start with a delay (contain a factor 1=z). The

readerwillrecognize

^ y(tjt 1)=W u (q)u(t)+W y (q)y(t) (14)

asthe Wienerlter, [18] for estimating (predicting)

y(t)fromu(s);y(s);st 1. Let

e

0

(t)=y(t) y(tjt^ 1)

Then(14)canberearrangedas

y(t)=G 0 (q)u(t)+H 0 (q)e 0 (t) (15a) with H 0 (z)=(I W y (z)) 1 (15b) G 0 (z)=H 0 (z)W u (z) (15c)

By the properties of the Wiener lter e

0

(t) will be

uncorrelatedwithy(s);u(s);st 1,i.e.

lim N!1 1 N N X t=1 e 0 (t)  y(t ) u(t )  =  0 0  ;8 1 (16) Sincee 0

(s)isconstructedfrom y(r);u(r);rs, this

alsoimpliesthat

lim N!1 1 N N X t=1 e 0 (t)e 0 (t )=0for 6=0 (17)

Thecorresponding limitfor  =0wedenote by

0 . Let (t)=  u(t) e 0 (t) 

  (z)=   u (z)  ue (z)  eu (z)  0  (18) where ue

(z)willbeananti-causalfunction,in view

of(16).

Remark. Notethat

ue

willnormallynotbezero,

evenif thereis nofeedbackin thedata. An explicit

example of a non-linear, causal, feedback-free

rela-tionship betweenu andy that still givesanon-zero

(butnon-causal)correlationbetweenuandeisgiven

inExample 1of[19]. .

The point of this discussion is of coursethat any

quasistationaryinput-outputdataset

fz(t);t=0;1;:::g

canbeseenasbeingproducedby(15a),withasignal

e

0

whichhasaconstantspectrum(\whitenoise")and

suchthat e

0

(t) is uncorrelatedwithpast u(s);s <t

(i.e. (16)holds.) Statisticalindependence betweene

anduandamongewillgenerallynothold. Anyway,

wehavenotintroducedanystochasticframeworkfor

thedata.

This means that considering just second order

properties (i.e. the spectra) of the signals y and u,

we cannot disprove that they have been generated by

(15a). In other words, the system (15a) is a

sec-ondorder equivalent ofthe systemthatgenerated

y fromu.

Now,itmustimmediatelybesaidthatG

0 andH

0

will in general depend on the input spectrum 

u ,

sothatthesecondorder equivalentobtainedforone

inputmaybeuselessto describethetruesystemfor

anotherinput.

4 A Characterization of the

Limit Model

We shall in this section develop some resultsabout

limits of estimated LTI-models based on data from

arbitrary systems. Thetheory will actuallybe

self-contained and it will not rely upon the traditional

convergence results for identied models, given e.g.

in[16].

Theresultisasfollows

Theorem 4.1 Consider the input-output data

fu(t);y(t);t = 1;2;:::g. Assume that the data are

quasistationary and that W

u

and W

y

given by (12)

{ (13) are well dened and stable. Consider the

LTI model structure (1) and let the estimate ^



N be

denedby (5a). Then

lim N!1 ^  N =argmin  Z   1 jH  j 2  (G 0 G  ) (H 0 H  )     u  ue  eu  0  (G 0 G  ) (H 0 H  )  d! (19) Here G 0 and H 0

are dened from u and y by (12)

-(15b), the argumente i!

of all the transferfunction

hasbeen omittedasin(2),andoverbardenotes

com-plexconjugation.

Notethatthisisexactlythesameresultthatholds

w.p.1in caseitisassumed that(15a) hasgenerated

the data with e

0

being a sequence on independent

randomvariancewithzeromeanvaluesandvariance

= 

0

. This is the basic, "traditional" convergence

result,see e.g. [16]. This meansthat alltraditional

analysisoflimitingestimatesinopenandclosedloop

can be directly applied to the general, non-linear,

non-stochastic case dealt with here, since that just

amountsto an analysis of theintegral in (19). See,

forexample,[20].

Toprovethistheoremwerstestablisharesultof

independentinterest:

4.2 Transformation ofSpectra by

Lin-ear Systems

Theorem 4.2 Letfw(t)gbeadeterministic,

quasis-tationarysignalwithspectrum

w

(!)andletG(q)be

astable lter. Let

z(t)= 

s(t)

w(t) 

isalso quasistationary withspectrum

 z (!)=  G(e i! ) w (!)G T (e i! ) G(e i! ) w (!)  w (!)G T (e i! )  w (!) 

The proof of this theorem is given in Appendix A.

Wemaynotethattheresultsstillparallelthetheory

of stationary stochastic processes. The expressions

fortransformingspectraareentirelyanalogous.

Forfamilies of linearlterswehavethe following

results.

Theorem 4.3 Let fG



(q); 2 Dg be a uniformly

stablefamilyoflinearlters(see(11))andletfw(t)g

beaquasistationary sequence. Let

s  (t)=G  (q)w(t) R s (;)= lim N!1 N X t=1 s  (t)s T  (t )

Then, forall 

sup 2D jj 1 N N X t=1 s  (t)s T  (t ) R s (;)jj!0asN !1 Proof

Weonlyhavetoestablishthat theconvergencein

(48)(in theappendix) to zerois uniform in  2D.

In the rst step all the g(k) terms carry an index

:g  (k). Interpreting g(k)=sup 2D jg  (k)j

(48)willof coursestill hold. SincethefamilyG

 (g) isuniformlystable 1 X k =0 g(k)<1

andthis wastheonly propertyoffg(k)g usedto

es-tablish that (48)tends to zero. This completes the

proof.



Thepredictionerrorsaccordingtothemodel(1)are

"  =H 1  (y G  u) (21)

where we havesuppressed all arguments. The

esti-mateisdeterminedbyminimization of

^  N =argmin N X t=1 " 2  (22)

Studying the second order properties of "



, we can

replaceywithitssecondorderequivalentdescription

(15a). Insertingthat expressionfory in(21)gives

"  =H 1  (G 0 u G  u+H 0 e 0 ) =H 1  [(G 0 G  )u+(H 0 H  )e 0 ]+e 0 M =v  (t)+e 0 (t) (23) AccordingtoTheorem 4.2",v  ande 0 are

quasista-tionarysignals,andaccordingtoTheorem4.3

N X t=1 " 2  (t)!  V()+ 0 (24) uniformlyin 2DasN !1 (25) where  V()= lim N!1 1 N N X t=1 v 2  (t) (26)

wherewealsousedthelimitsin(16)and(17). With

thenotation of(7b) 

V()=R

v



(0), sofrom the

in-verse Fourier transform (or Parseval's relationship)

wehavethat  V()= Z    v (e i! )d! (27) where v isthespectrumofv  ,whichaccordingto

(23)andTheorem 4.2isgivenby

1 jH  j 2  (G 0 G  ) (H 0 H  )   (28)   u  ue  eu  0  (  G 0  G  ) (  H 0  H  )  (29)

Any estimated model will be an imperfect

descrip-tionofthesystem. ThetermModelErrorModelwas

coined in [4] to denote any way to characterize the

errorsassociatedwiththemodel. These wayswillof

coursethemselvesbeimperfect,buttheymaybe

ad-equatetodescribetheamountofcautionthatshould

be exercisedwhen the nominal model is used. The

basic model error model could simply be described

byaparallelblocktothenominalmodelisshownin

Figure7.

How do we gaininformation aboutthe model

er-ror model? Well, all information is in the

mea-sureddata, possiblyin conjunction withsome

data-independent prior knowledge. Since the nominal

model has squeezed out most{ orpart{ of the

in-formationinthedata,themodelerrormodelwill

de-scribe therelationshipbetweentheinput uand the

outputerrorv(t)=y(t) G(q; ^



N

)u(t)orthe

resid-uals "(t) = "(t; ^



N

). This is also illustrated in

Fig-ure7. Consequently,developingamodelerrormodel

amounts to some kind of residual analysis. This is

a standard topic in regression theory, see e.g. [21],

and the analysis of correlation between past inputs

and residuals,depicted, e.g. in Figure 3isthemost

commonexampleofsuchanalysis.

v y G mem u u G nom

Figure 7: Thenominal model G

nom

and themodel

errormodelG

mem .

Building linearmodelerrormodelsisthus justan

alternative wayof phrasing the resultof such

stan-dard (second order) residual analysis. See [4].

Ex-plicitlinearmodelerrormodelswillconsequently

de-scribethebiasdistributionofthenominalmodel,but

willhavenoinformationaboutpossibleerrorsdueto

Thenominalmodelplusthelinearmodelerrormodel

willjust describetheLTI-equivalent,dened in

Sec-tion3.

Itisthereforeofmoreinteresttodiscusserror

mod-els that are nonlinear and/or time-varying. A brief

discussionof this isgiven in [22]. Now,the purpose

of an error model is not to complement the

nomi-nalmodelwithdetailedstructuralinformation. That

shouldratherbedoneaspartofthenominalmodel.

Instead,thepurposeoftheerrormodelistocapture

the reliability of the nominal model, so that proper

robustnessin thecontroldesigncanbeassured.

Thismeansthatweshallworkwithamodelerror

model depicted in Figure 8. We shall only be

con-cerned about the gain of the block ~g

mem

. Written

outasequationswehave

" v ~ g mem u F u W 1 W 2

Figure8: Themodelerrormodelwithlinear

weight-ingfunctions "(t)=W 1 2 (q)v(t) (30a) u F (t)=W 1 (q)u(t) (30b) "(t)=g~ mem (u t 1 F ) (30c) k"kku F k+ (30d)

Somecommentsareinorder:

 Theroleof theweightingfunctionsW

1 andW

2

is to give adequatefreedom for the control

de-sign. Estimating just the gain of the middle

block could be an obtuse instrument, and the

linearweightswillproveuseful.

 Thenorms in (30d)are to beinterpreted in L

2

sense.With=0,thenumberisconsequently

theH

1

gainofthesystemg~

mem .

1. Toallowforexternalsignalstoenterthe

er-rormodel,asdepictedinFigure9(would

thenbethenormofw)

2. To allow for possible very large gains for

smallamplitudesignals,whichmaynotbe

harmful for "practical stability". This is

further elaborated in [23]. Fordiscussions

of such ano-set term in connection with

stabilityseealso[24]and[25].

w " v ~ g mem u F u W 1 W 2

Figure 9: Themodelerror model with additive

dis-turbance

6 Estimating the Gain of a

Sys-tem

Weare nowfacedwith anessentialproblem: Given

the sequencesu

F

and",how toestimate  andin

(30d)?

There isapparentlynotanextensiveliteratureon

this problem. Some "identication for robust

con-trol" articles relate to the gain estimation, like [6],

[7], [26],[27], [28] and [29]. These mostly dealwith

thegainofaLTIoranLTVerrormodel,though.

It is not the purpose of this section to launch a

recommendedmethod forgainestimationofgeneral

model errormodels. Weshall insteadpointto some

possibilities, that indicate that the problem is not

infeasible.

A ratherobviouspossibilityis to explicitlyestimate

themodelin(30c)andthencomputethegainofthe

estimatedmodel: "(t)=g~ mem (u t 1 F )+w(t)

Use yourfavoritenon-linear black box model

struc-turefor~g

mem

, such asanArticialNeuralNetwork,

Local Linear Models, Piece-wise linearmodels, etc.

(cfChapter5in[16]). Thendetermine andfrom

theestimatedmodelandthesizeofw.

Asan alternative,ifjust the gainis of interest,it

may be simpler to directly estimate a "ceiling" for

thesurfacethatg~

mem

denes. SeealsoSection6.3.

6.2 Estimating the Gain Directly

From Data

Itistemptingtocircumventthelaboriousprocessof

estimating ageneralnonlinearblack-box model and

thencomputeitsgain,bydirectlyestimatingthegain

from the data. Forexample, if a local, radialbasis

neuralnetworkisusedtoestimatethesurface~g

mem ,

the"peaks"ofthissurfacearecreatedbylargevalues

ofobserved"(t). (SeeSection 6.3 formoreintuition

aboutthis\surface.") Thehighestgainpointsofthe

surface are created by observations where the ratio

j"jtokukislarge. Thisleadstothefollowingsimple

method:

 Assume that it is known that most of the

in- uenceon"(t) from pastu

F

(s);st,linear or

not,lastsfordsamples. Simpletransient

exper-iments,orbasicpriorknowledgecangiveinsight

intothis. Form

'(t)=  u F (t 1) u F (t 2) ::: u F (t d)  (31a) andnd =max t j"(t)j k'(t)k (31b) Here k
k 2

is the usual 2-norm. Now,  is the

seen in the data. To move to a corresponding

norm for " anatural upper bound on the gain

wouldbe

^

 = p

d
 (31c)

Thereasonwhy ^

 isanupperbound,isthatitdoes

not follow that dsuch large valuesof " canbe

pro-duced in a sequence. Now this is averysimple

al-gorithm,that doesnothaveanyprovisionsfor

deal-ing with noise or o-sets. A more general version

wouldbetohaveanintelligentwayofndinga

noise-permissive upper-bounding line when regressing j"j

onk'k. Herewejustletthatlinegothroughthe

ori-gin(=0)anddidnotallowanyobservationsabove

theline.

Anyway,letustesthowthisestimatorworks.

Example 6.1: Estimating Gains for

Time-Varying, Non-Linear, Noise Corrupted

Sys-tems

We create a time-varying, non-linear, noise corrupted

systemasfollows:

 Create two random, linear third order system:

m1=idpoly(fstab([1,randn(1,3) *2], ...

[0,randn(1,3)*3]) andsimilarlyform2.

 Create an input signal u as a white noise normal

signal with 1000 samples and low pass lter it by

1=(q 0:8)

 Let u pass through a static, discontinuous

non-linearitytoformu 1 : u 1 = ( 5u ifjuj2 u else

 Formatimevaryinglinearsystemfromm1andm2by

letting itsparameters varyas acosine with period

200samplesbetweenthoseofm1andm2. Theoutput

whensimulatedwithu1 iscalledy1.

 Introduceanoutputdead-zonesothat

y(t)= ( 0 ifjy 1 (t)j5 y 1 (t) else

 add rectangular distrubuted noise to y so that the

signal-to-noiseratiobecomes10(amplitude-wise)

0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

80

100

120

140

160

180

200

Figure10: Evaluationofthegainestimator(31). The

plot shows the result for 200 simulated systems as

describedinthetext. Eachdotcorrespondstoa

sys-tem. Its y-coordinate is the estimated gain and its

x-coordinateisthetruegain.

Thetheoreticalgainofthis,non-linear,time-varying

sys-temis 5times thelargestmagnitudethatthe frequency

functionsofm1andm2everassume. Twohundredsystems

ofthiskindweresimulated. Onlysystemsm1andm2with

impulseresponsesolutiontime(to5%)lessthan20

sam-pleswereaccepted. Thereasonisthatsystemswithlong

impulse responses probably require modied techniques

forgainestimation(seeSection6.3).

Figure10showsthegainestimatefromalgorithm(31)

versusthetruegain. Therootmeansquaredeviationof

themeasure

EstimatedGain

TrueGain

1 (32)

is34%,whichcouldbeperceivedasasurprisinglygood

result.

6.3 General Gain Estimates: A

Dis-claimer

Itisinstructivetovisualizethegainestimation

prob-lem asfollows: ConsiderR d+1

. Letthe \Floor"R d

usview"asrisingperpendiculartothis oor. A non-linearmodel~g mem asin(30c)thenisahypersurface overR d

. Estimatingthe gainis amatterof nding

thehighestelevations ofthis surfaceasviewed from

theorigin.

Now,R d

is aprettybigand \empty" space.

Sup-poseweuse d=20 asin theexample. Considerthe

unit cube ju

F

(t)j 1and usea gridof neness 0.2

to distinguishbetweenvaluesof u

F

, whichis rather

crude. Then the unit cube will contain 10 20

cells.

Even with quite a respectable number of

observa-tions, like N = 10 4

, at most a portion of 10 16

of

the cells will be populated with observations. The

surface mentioned above will therefore have an

ex-tremely thin support of observations. Finding, and

estimatingtheangletothepeaksofthissurface

con-sequently will be a tricky problem. Practically

re-gardless of the number of observations made, most

parts ofthe spacehavenotbeencovered,and

with-outpriorinformationitisimpossibletosaywhatthe

gainwouldhavebeenatthoseparts.

ItisinthelightofthisthattheresultsofFigure10

couldbeconsideredas \surprisinglygood".

Now,thelongertheeectofaninputsamplelasts,

the morediÆcult will the gainestimation be, since

theprobabilitywewillhitthe\worstcase"input

se-quence becomes less. This was the reason that we

onlystudiedsystemswithsolutiontimelessthan20

samples,whichanywayisareasonablylongresponse

time. Systemswith longer lasting responses will

re-quiremodiedestimationtechniques.

What can be done about this lack of support of

observationsinR 20

? Well,essentiallynothing. Some

possibilitiesto dealwiththeproblemcouldbe:

1. Obtainmoremeasurements: Willnothelpmuch,

since R 20

would require a totally unrealistic

amountofdatato becovered.

2. Assume that the surfaceis\verysmooth", and

that the collected data exhibit the behavior of

thesystem,that weare likelyto encounter also

later. This is really the alibi behind algorithm

(31).

3. Assume that the surface is a hyperplane, i.e.

mem

that theactualmodel surfacecanbeeectively

over-boundedbysuchahyperplane.

4. Assume that the surface can be well

approxi-matedbyaradialbasisneuralnetwork. This is

essentiallythesameas2.

5. Assume that the surface can be well

approxi-matedwitharidgetypeneuralnetwork,suchas

the traditional sigmoidal networks. This, in a

sense, is a combination of 2 and the idea that

youcanextrapolatealonghyperplanes.

It is obvious, in thelight of this, that noprocedure

for estimating the gain can come with any quality

guarantees, unless someveryreliableprior

informa-tionisavailableabouttheshapeofthe surface. For

alinearerrormodel,itwouldbepossibletodescribe

thedistributionoftheestimateprovidedby(31),but

in thegeneralcasesuchanalyticalresultscannotbe

derived.

Estimatingthegaininthegeneralcasewillthusbe

subject to verication in the particular applications

of interest, just as the constructionof general

non-linearblack-boxmodels.

7 An LTI Model with a

Gen-eral Model Error Model as

an Equivalent Uncertain LTI

Model

7.1 Linear Model Errors

Onceamodelwithitsmodelerroruncertaintyis

de-livered,thequestionishowtodesignacontrollerthat

will stabilizethe systemrobustly. Bythis wewould

mean that the chosencontroller should stabilize all

models in the \region"denedby the nominalmodel

andthe model errormodel.

Incase we haveused a linearmodel errormodel,

thisregioniseasilydepictedinthefrequencydomain.

Itwilllooklikeastripin theBode,orNyquistplot,

i.e. G2G=fGj jG(e i! ) G nom (e i! )j<
(!)g (33)

for such a set of models is well known: Choose a

regulatorK,suchthatthecomplementarysensitivity

function T = G nom K 1+G nom K (34)

islessthantheinverserelativemodelerrorbound:

jT(e i! )j< jG nom (e i! )j
(!) ; 8! (35) H 1

techniquescanbeusedtodetermineifsuchaK

exists,forgivenG

nom

and
. See,e.g. [30].

7.2 Frequency Weighted Non-linear

Model Error Model

w " v ~ g mem G nom u F u y K W 1 W 2

Figure 11: Block diagramof thefeedback loopwith

modelerror

Theerrormodel(30)correspondsto aclosedloop

block diagram as in Figure 11. This can be

rear-rangedtobeseenasfeedbackbetweenthenon-linear

partof theerrormodelg~

mem and KW 1 W 2 1+KG nom

(keepinginmindthatweonlyconsiderSISOmodels

here). Suppose that the gainof the non-linear part

issubjectto

k"kku

F

k+ (36)

eects of the non-linearity and of the additive

dis-turbance w. The small gain theorem tells us that

stabilityisassuredif     W 1 (e i! )W 2 (e i! )K(e i! ) 1+K(e i! )G nom (e i! )     <1 8! (37)

Comparingwith(35)werealizethatwejustcan

con-sider the set of possible system descriptions to be

linearand givenby

G2G=fGj jG(e i! ) G nom (e i! )j<::: (38) <W 1 (e i! )W 2 (e i! )g

By stabilizing any linear model in this set, i.e.,

achieving (35) for
=W

1 W

2

, wehave also made

the linear control design robust against non-linear

model errors ofthe type (30).

We can also go beyond stability robustness and

consider sensitivity to disturbances. It follows, see

[23],that theoutputnormis bounded by

kykkSW 2 k  1 kG w k ; G w =T W 1 W 2 G nom (39)

whereS isthesensitivity andT thecomplementary

sensitivity of the nominal design. Again, standard

lineartechniquestellushowtodesignthepairSand

T fromW 1 ;W 2 ;; andG nom

sothatthesensitivity

expressedby(39)isacceptable.

7.3 An Equivalent Uncertain Linear

ModeltobeDeliveredtothe User

From the discussion above it follows that if the

LTIidenticationprocessestimatesanominalmodel

G

nom

andweselecttheweightingfunctions W

1 and

W

2

and thenestimate thegain of theblockg~

mem

we can deliver an LTI uncertainty model consisting

of G

nom

and the band G dened by (38). If robust

linear control design is applied to this uncertainty

model, LTI regulators will beproduced that are

ro-bustalsotonon-linear,time-varyingmodelerrorsup

to the size determined by the gain estimator. This

extends in aquite naturalway LTI-identication +

LTIcontroldesigntogeneralsystemsthatcanbewell

It may be quiteimportantto correctly usethe

free-dom oeredby theweightsW

1

and W

2

. As will be

seeninthenextsection,dierentweightscanproduce

quite dierent LTI uncertainty models. The choice

ofW isaninterplaybetweenshapingtheuncertainty

regions to what suits the control design, and

creat-ingdescriptionsthat leavetheunexplained (\") as

smallaspossible.

Somenaturalchoices are

 W 1 =G nom . This makesu F =y,^ themodel's

simulatedoutput. It is naturalto comparethe

modelerrorwiththesimulatedoutput,sincethis

directlyrelatestothepercentageoftheoutput's

variationthat isexplainedbythemodel. Italso

leadstoaquanticationoftherelativemodel

er-ror,whichnaturallyarisesinrobustnesscriteria

(see e.g. G

w

in (39) which contains the ratio

W 1 =G nom ).  W 2 = H nom

, the nominal noise model. This

makes"equaltothemodelresiduals,whichgives

anoutputthetheunknownblockwiththe

small-est possible variance. This should leadthe the

smallest,but theshapeoftheuncertainty

re-gionmayperhaps be unsuitable for control

de-sign.

Thereareofcoursemanyotherpossiblechoices. One

should however avoid weighs with long impulse

re-sponses, since this may make the gain estimation

moretricky.

8 Some Numerical

Experimen-tation

Letus dosomeexperimentsto seehowtheoutlined

works out. We rst test a time-varying, nonlinear

system:

Example 8.1: Estimating LTI models for

Non-Linear, Time-Varying Systems

Considerasystemthatistime-varyingbetweenthetwo

y(t) 2y(t 1)+1:45y(t 2) 0:35y(t 3)

=u(t 1)+0:5u(t 2)+0:2u(t 3)

and

y(t) 1:93y(t 1)+1:43y(t 2) 0:41y(t 3)

=1:05u(t 1)+0:41u(t 2)+0:18u(t 3)

Itisalsosubjecttoaninputstaticnon-linearity,so that

inputswith anamplitude less than0.8 is multipliedby

1.2,aswellasanoutputdead-zoneoflength1.Theinput

is whiteGaussian noisewith unitvariance. Athird

or-derLTImodelwasestimatedfromthedata. Thispasses

the traditional model validations tests well. Figure 12

shows the nominal estimatedmodel and the equivalent

uncertainLTI-models,as described inthe previous

sec-tion,withthegainestimatedusing(31).

Finally, wereturnto Example1.1.

Example 8.2: Rotation of a Rigid Body,

Cont'd

FromthedataofexperimentsAandB(seeFigures1

-2)nominalthirdorderLTImodelswereestimatedas

de-scribedinExample1.1. Error modelswereestimatedas

in(30)and(31)forsomedierentW1andW2. Figures13

and14showtheamplitudebodeplotsoftheresulting

er-rormodels.TheseshouldbecomparedwithFigures5and

6. We seethat the essentially linearcase of experiment

Aiscorrectlyidentiedassuch,whilethenon-linearcase

ofexperimentBgivesanerrormodelthatclearly shows

thatareliablelinearapproximationisnotfeasible.

9 Conclusions

Inthiscontributionfourfactshavebeenpointedout:

 Undergeneralconditionswecanexplicitly

spec-ifyinwhichwayanestimatedLTImodel

approx-imatesageneralsystem. Itisessentiallyonly

re-quiredthatthesystemproducesquasistationary

10

−2

10

−1

10

0

10

1

10

−1

10

0

10

1

10

2

Frequency (rad/s)

Amplitude

From u1 to y1

10

−2

10

−1

10

0

10

1

10

−1

10

0

10

1

10

2

Frequency (rad/s)

Amplitude

From u1 to y1

10

−2

10

−1

10

0

10

1

10

−1

10

0

10

1

10

2

Frequency (rad/s)

Amplitude

From u1 to y1

10

−2

10

−1

10

0

10

1

10

−1

10

0

10

1

10

2

Frequency (rad/s)

Amplitude

From u1 to y1

10

−2

10

−1

10

0

10

1

10

−1

10

0

10

1

10

2

Frequency (rad/s)

Amplitude

From u1 to y1

Figure 12: Results from the experiment described

in Example 6. The amplitude bode plots show as

alightshadedregiontheerrormodelsconstructedas

in Section4. Thedarkshadedregionis thenominal

estimated LTI model along with an uncertainty

re-gioncorrespondingto1standarddeviation. Thefour

thin linesarethefrequencyfunctions ofthetwo

lin-earsystems,eachmultipliedby1andby1.2(Recall

thatthereisastaticnon-linearitywithgainbetween

1and1.2.) Theplotscorrespondtodierent

weight-ingltersW

1 andW

2

. Fromaboveandleftto right:

(1)W 1 =W 2 =1, (2)W 1 =G nom ,W 2 =1, (3)W 1 =G nom ,W 2 =H nom

(nominalnoisemodel).

(4)W 1 =1=(q+0:3);W 2 =1, (5)W 1 =1=(q 0:95);W 2 =1.

10

0

10

1

10

2

10

−5

10

−4

10

−3

10

−2

Frequency (rad/s)

Amplitude

From u1 to y1

10

0

10

1

10

2

10

−5

10

−4

10

−3

10

−2

Frequency (rad/s)

Amplitude

From u1 to y1

Figure 13: The resultinguncertainty model for

Ex-perimentAinExample1. Left: Relativemodelerror

(i.eW 1 =G nom )withW 2 =H nom . Right: Relative

modelerrorwithW

2 =1.

10

0

10

1

10

2

10

−5

10

−4

10

−3

10

−2

Frequency (rad/s)

Amplitude

From u1 to y1

10

0

10

1

10

2

10

−5

10

−4

10

−3

10

−2

Frequency (rad/s)

Amplitude

From u1 to y1

Figure 14: The resultinguncertainty model for

Ex-perimentBinExample1. Left: Relativemodelerror

(i.eW 1 =G nom )withW 2 =H nom . Right: Relative

modelerrorwithW

2 =1.

estimating the size of thedistance betweenthe

truesystemandtheLTI-approximation

 Wehaveshownhowtheresultingmodelcanbe

seenasanLTI-modelwithanuncertaintyregion,

muchinthesamespiritasthetraditionalmodel

withstatisticalcondenceintervals.

 LTI robustcontrol design forthe familyof LTI

models delivered by this process will give

reg-ulators that are robustalso to model errors

re-sultingfromthepossiblynonlinear,time-varying

truesystem

An artifactofthestandardLTI identication

ma-chinery is that it produces a nominal model with a

condence interval that tends to zero as the

num-berof observeddata growsto innity. This isreally

anundesiredfeature,since,realistically,thereareno

trueLTIsystemsin thereal world.

An attractiveaspect ofthe outlinedwayof

deliv-ering uncertainLTI models is that it resembles the

classicalapproach,withtheimportantexceptionthat

theuncertaintyregionswilltypicallynottendtozero

asmoreandmoredatabecomeavailable. Therewill

be some \remaining uncertainty", which should be

thoughtofasahealthysign.

Now, the outlined process also will need several

enhancements:

 More eective gain estimators are required.

Thereshould beagoodpotentialforsuch a

de-velopment. Thefundamental limitation is that

youcanonlybasetheestimateonwhatyouhave

seenandtypicallytheobservationsarebutatiny

fraction of the actual response surface. This is

morepronouncedifthe response timeto an

in-putchangeislong. Theneedtodealwithworse

signal-to-noiseratiosthanthatinFigure10calls

fortechniquesthatallowcertainobservationsbe

outsideaboundingconeorabounding\ceiling"

of the response surface. For a time-invariant

system this should be quite feasible, but for a

time-varyingsystemthedistinctionbetween

sig-nalandnoiseisnottrivial.

servative. Thisisnotjustaconsequenceofpoor

gainestimates, but another reason isthat

hav-ingjust againmeasurewill notrevealmuch of

thestructureoftheuncertainty. Putdierently,

thesmallgaintheorem isquiteconservative. It

wasillustratedinFigure12howtheuncertainty

regionsmaydependonthechosenweightsinan

essentialway. Amoregeneralerrormodelwould

betoestimatethegainforablock

u F =W 1 u+W 12 v to "=W 21 u+W 1 2 v (40)

This corresponds to an error model as in

Fig-ure15,whichiswellpreparedforLTIcontrol

de-W " ~ g mem u F u v

Figure15: Amoregeneralmodelerrormodel. The4

transferfunctionsin thelinearblock W arerational

combinations of the functions W

1 ;W 2 ;W 12 ;W 21 in (40).

sign,usinge.g. H

1

techniques. Thecasein

Fig-ure9clearlyisthespecialcaseW

12 =W

21 =0.

Thetwoextraweightingfunctionswillgivemore

freedom to customize theerror description. At

the same time, the resulting LTI uncertainty

model (consisting of G

nom

, the four transfer

functionsin W and thegainestimate) isnow

notsimplyaband aroundthe Nyquistcurveof

G

nom .

 Athirdlineofthoughttopursue,istomovefrom

symmetric descriptions, using e.g. IQC's, [31],

[32]. Whilethegainestimatein(31)amountsto

ndingthescalar in expressionslike

Z  "(t) u F (t)   1 0 0  2  "(t) u F (t)  dt 0 8u F ;" (41)

whichalsocanbewrittenintermsoftheFourier

transformsofthesignals. Themoregeneralcase

(40)correspondsto Z  V( i!) U( i!)   "  2 jW 12 j 2 jW 2 j 2  2 W 1 W 12 W 21 W 2 1  2 W 1 W 12 W 21 W 1 2  2 jW 1 j 2 jW 21 j 2 #   V(i!) U(i!)  d!0 8u;v (42)

The IQC approach would be to nd a matrix

(!)suchthat Z  V( i!) U( i!)  (!)  V(i!) U(i!)  dt 0 8u;v (43)

The kinship with the gain estimation is clear

from (43), (42). In this case,the delivered LTI

uncertainty model would be fG

nom

, g which

maycontainmorestructural informationabout

thecharacter oftheuncertainty, relatedto

pas-sivity properties. Control design basedonsuch

anuncertaintymodelisdiscussed, e.g. in[32].

Acknowledgments Thisworkhasbeensupported

bytheSwedishResearchCouncil(Vetenskapsradet).

DiscussionswithTorkelGladandAndersHelmersson

(who together also suggested Example 1.1), as well

aswithWolfgangReineltandAlexanderNazinhave

beenveryhelpful inthepreparationofthisarticle.

4.2

Proof

Firstassumethatw(s)=0fors0andconsider

R N s ()= 1 N N X t=1 s(t)s T (t ) = 1 N N X t=1 t X k =0 t  X `=0 g(k)w(t k)w T (t  `)g T (`) (44)

With the conventionthat w(s)=0if s62 [0;N] we

canwrite R N s ()= N X k =0 N X `=0 g(k)  1 N N X t=1 w(t k)w T (t  `)g T (`) (45) Let R N w ()= 1 N N X t=1 w(t)w T (t ) WeseethatR N w

(+` k)andtheinnersumin(45)

dierbyat mostmax(k;j+`j) summands, each of

which arebounded byC accordingto (7a). Thus

jR N w (+` k) 1 N N X w(t k)w T (t  `)j C max(k;j+`j) N  C N (k+j+`j) (46) Letusdene R s ()= 1 X k =0 1 X `=0 g(k)R w ( +` k)g T (`) (47)

jR s () R N s ()j  X 0 X 0 jg(k)jjg(`)jjR w (+` k)j + N X k =0 N X `=0 jg(k)jjg(`)j jR w ( +` k) R N w ( +` k)j + C N N X k =0 kjg(k)j
N X `=0 jg(`)j + C N N X `=0 j+`jjg(`)j
N X k =0 jg(k)j: (48)

Here,therstsumisoverthecomplementaryindices

ofthesecondonei.e. k>Nand/or`>N. Thisrst

sumtendsto zeroas N!1sincejR

w

()jC and

G(q) is stable. It follows from the stability of G(q)

that 1 N N X k =0 kjg(k)j!0asN !1 (49)

Hencethelasttwosumsof(48)tendtozeroasN !

1. Considernowthesecond sumof(48). Selectan

arbitrary">0andchooseN =N

" suchthat 1 X k =N"+1 jg(k)j<"=[C
C 1 ] (50) where C 1 = 1 X k =0 jg(k)j

ThisispossiblesinceGisstable. ThenselectN 0 " such that max 1<`<N " 1<k <N " jR w (+` k) R N w (+` k)j<"=C 2 1 forN >N 0 "

. Thisispossiblesince

R N w ()!R w ()asN !1 (51)

ber (which depends on ") of R

w

(s):s are involved

(nouniformconvergenceof(51)isnecessary). Then

for N >N 0

"

we havethat thesecond sum of(48) is

bounded by N X k =0 N " X `=0 jg(k)jjg(`)j
" C 2 1 + 1 X k =N " +1 1 X `=0 jg(k)jjg(`)j
2C + 1 X k =0 1 X `=N"+1 jg(k)jjg(`)j
2C

which is lessthan 5" accordingto (50). Hence also

thesecondsumof(48)tendstozeroasN !1,and

wehaveprovedthatthelimitof(48)iszero,andthat

hences(t)is quasistationary.

Theproofthat lim(1=N) P

N

t=1

s(t)w(t ) exists

isanalogousandsimpler.

For s (!)wenowndthat  s (!)= 1 X = 1 1 X k =0 1 X `=0 g(k)R w (+` k)g T (`) ! e i! = 1 X = 1 1 X k =0 g(k)e ik !  1 X `=0 R w ( `+k)e i(+` k )! g T (`)e i`! =[ `+k=s] = 1 X k =0 g(k)e ik !
1 X s= 1 R w (s)e is!
1 X `=0 g T (`)e i`! =G(e i! ) w (!)G T (e i! )

Hencetheupperleft cornerof

z

(!)isproven. The

odiagonaltermsareanalogousand simpler.



References

[1] B. Wahlberg and L. Ljung, \Design variables

forbiasdistributionintransferfunction

estima-tion,"IEEETrans.AutomaticControl,vol.

modelquality," Automatica,vol.31,no.12,pp.

1771{1797,1995.

[3] G.C. Goodwin, M. Gevers, and B. Ninness,

\Quantifying the error in estimated transfer

functionswithapplicationtomodelorder

selec-tion," IEEE Trans.AutomaticControl,vol.37,

no.7,pp.913{929,1992.

[4] L. Ljung, \Model validation and model error

models,"inThe 

AstromSymposiumonControl,

B.WittenmarkandA.Rantzer,Eds.,pp.15{42.

Studentlitteratur,Lund,Sweden,August1999.

[5] L.LjungandL.Guo,\Theroleofmodel

valida-tionforassessingthesizeof theunmodeled

dy-namics," IEEE Trans. Automatic Control, vol.

AC-42,pp.1230{1240,1997.

[6] R. Smith and J. C. Doyle, \Model validation:

aconnectionbetweenrobustcontroland

identi-cation," IEEE Trans.AutomaticControl, vol.

AC-37,pp.942{952,1992.

[7] R.SmithandG.E.Dullerud, \Continuous-time

controlmodelvalidationusingnite

experimen-taldata," IEEE Trans.AutomaticControl, vol.

AC-41,pp.1094{1105,1996.

[8] R. Kosut, M. K. Lau, and S. P. Boyd,

\Set-membershipidenticationofsystemswith

para-metric and nonparametric uncertainty," IEEE

Trans.AutomaticControl, vol.AC-37,pp.929{

941,1992.

[9] R.G. Hakvoort and P.M. van den Hof,

\Iden-ticationofprobabilisticsystemuncretianty

re-gionsbyexplicitevaluationofbiasandvariance

errors," IEEE Trans.Autom. Control, vol.

AC-42,no.11,pp.1516{1528,Nov1997.

[10] MichelGevers, \Towardsajointdesignof

iden-tication and control?," in Essays on control:

Perspectives in the theory and its applications,

H L Trentelman and J C Willems, Eds., ECC

'93Groningen,1993.

measureofrobuststabilityforanidentiedsetof

parameterizedtransferfunctions," IEEETrans.

Autom. Control, vol. AC-45, no. 11, pp. 2124{

2145,Nov2000.

[12] HakanHjalmarsson,AspectsinIncomplete

Mod-eling in System Identication, Ph.D. thesis,

LinkopingUniversity,Linkoping,Sweden,1993,

Linkoping Studies in Science and Technology,

DissertationsNo.298.

[13] P. Van Overschee and B. DeMoor, Subspace

Identication of Linear Systems: Theory,

Im-plementation, Applications, Kluwer Academic

Publishers,1996.

[14] M.Verhaegen, \Identication of the

determin-isticpartofMIMO statespacemodels,givenin

innovationsform frominput-outputdata,"

Au-tomatica,vol.30,no.1,pp.61{74,January1994.

[15] W.E. Larimore, \Canonicalvariateanalysisin

identication,lteringandadaptivecontrol,"in

Proc. 29th IEEE Conference on Decision and

Control, Honolulu,Hawaii,December1990, pp.

596{604.

[16] L.Ljung, SystemIdentication-Theory forthe

User, Prentice-Hall, UpperSaddle River, N.J.,

2ndedition,1999.

[17] L.Ljung,\Anonprobabilisticframeworkfor

sig-nalspectra," in Proc. 24th IEEE Conf. on

de-cision andControl, FortLaudredale, Fla,1985,

pp.1056{1060.

[18] Norbert Wiener, The Extrapolation,

Interpo-lation and Smoothing of Stationary Time

Se-rieswith Engineering Applications, Wiley,New

York,1949.

[19] Urban Forssell and Lennart Ljung, \A

pro-jection method for closed-loop identication,"

IEEE Transactions on Automatic Control, vol.

AC-45,no.11,pp.2101{2106,Nov2000.

[20] U. Forssell and L. Ljung, \Closed loop

identi-cation revisited," Automatica, vol. 35, no. 7,

Analysis,2nded., Wiley,NewYork,1981.

[22] Lennart Ljung, \Model error modeling and

control design," in Proc IFAC Symposium

SYSID2000.,SantaBarbara,CA,Jun2000,pp.

WeAM1{3.

[23] S.T. Glad, A. Helmersson,and L. Ljung,

\Ro-bustness guarantees for linar control designs

withanestimatednonlinearmodelerrormodel,"

in IEEE Conference on Decision and Control,

Orlando,2001, Submitted.

[24] V.A. Yakubovich, \Absolute stability of

non-linear control systems in critical cases.,"

Av-tomatika i Telemekhanika, vol. 24, no. 3, pp.

293{302,1963.

[25] M. Vidyasagar, Nonlinear Systems Analysis,

Prentice-Hall,EnglewoodClis,NJ,1993.

[26] R.A. Davis and K. Glover, \An application of

recentmodelvalidationtechniquesto ighttest

data," in Proc. of the Third European Control

ConferenceECC95,Rome,Sept, 1995,pp.1249

{1254.

[27] R. L. Kosut, \Uncertainty model

unfalsica-tion," in IFAC Symposium on System

Identi-cation, Santa Barbara,CA,2000.

[28] R.L.Kosut, \Iterativeadaptiverobustcontrol

viauncertainty modelunfalsication," in Proc.

13th IFAC Congress, San Francisco,July 1996,

International federation of Automatic Control,

pp.91{96.

[29] K.Poolla,P.Khargonekar,A.Tikku,J.Krause,

and K. Nagpal, \A time domain approach to

modelvalidation,"IEEETrans.Automatic

Con-trol, vol.AC-39,pp.951{959,1994.

[30] K. Zhou, J. C. Doyle, and K. Glover, Robust

andOptimalControl, Prentice-Hall,Englewood

Clis,NJ,1996.

[31] V.A. Yakubovich, \The method of matrix

in-equalitiesin thetheory ofstabilityofnonlinear

no.7,pp.1017{1029,1964.

[32] A.MegretskiandA. Rantzer, \Systemanalysis

viaintegralquadraticconstraints,"IEEETrans.

Autom. Control,vol.AC-42,no.6,pp.819{830,