16.05.2016

1

prin persistenţa semnalului de intrare. A doua condiţie fiind

îndeplinită numai dacă v(k) este zgomot alb, constituie o

limitare a metodei.

Variabilele instrumentale sunt introduse tocmai în ideea de

a obţine un estimator similar celor mai mici pătrate care să fie

consistent indiferent de natura perturbaţiei. Metoda de variabilă

instrumentală utilizează estimatorul

N

k

N

k

TVI kykz

Nkkz

N 1

1

1

)()(1

)()(1ˆ (5.76)

în care (k) are aceeaşi semnificaţie ca în cazul metodei CMMP ,

Tnbkukunakykykyk )(),....,1(),(),......,2(),1()(

iar z(k), de dimensiune nbnakz )(dim este un vector oarecare

ale cărui componente, numite variabile instrumentale, trebuie

alese în aşa fel încât estimatorul să fie consistent

Înlocuind )()()( kvkky T în (5.74), se obţine

N

k

N

k

T

N

k

N

k

T

N

k

TN

k

TVI

kvkzN

kkzN

kvkzN

kkzN

kkzN

kkzN

1

1

1

*

1

1

1

*

1

1

1

)()(1

)()(1

)()(1

)()(1

)()(1

)()(1ˆ

(5.77)

Din (5.77) rezultă condiţiile de consistenţă a estimatorului VI:

N

k

TT kkzN

kkzER1

0)()(1

)]()([

N

k

kvkzN

kvkzE1

0)()(1

)()( (5.79)

(5.78)

Variabilele instrumentale pot fi corelate cu intrările şi ieşirile, dar nu sunt

corelate cu perturbaţiile. Pentru a satisface cerinţele (5.78), (5.79) foarte

frecvent se aleg ca variabile de stare eşantioanele întârziate ale intrării. La

întârzieri mai mari condiţiile sunt mai bine satisfăcute .

16.05.2016

2

Ca şi estimatorul CMMP şi estimatorul VI este insensibil la o

transformare liniară. Această proprietate se poate utiliza la

construirea vectorilor de variabilă instrumentală. Într-adevăr, dacă în

(5.76) înlocuim z(k) cu Tz(k) , unde T este o matrice n x n de

transformare nesingulară, [28], [77] estimatorul devine

N

k

N

k

TVI kykTz

NkkTz

NT

1

1

1

)()(1

)()(1ˆ

VI

N

k

N

k

T

N

k

N

k

T

kykzN

kkzN

kykzN

TTkkzN

ˆ)()(1

)()(1

)()(1

)()(1

1

1

1

1

11

1

5.5.2. Alegerea variabilelor instrumentale de bază

Se consideră o formă generală a vectorului de variabilă instrumentală

T

nbkukunakkqKkz )(),....,1(),(),...,1()()( 1 (5.80)

unde (k) este obţinut prin filtrarea datelor de intrare

[28],[76] :

ndnd

ncnc

qdqdqD

qcqccqC

kuqD

qCk

....1)(

....)(

)()(

)()(

11

1

110

1

1

1

(5.81)

K(q-1) şi K-1(q-1) sunt filtre asimptotic stabile, iar polinoamele

C(q-1) şi D(q-1) au zerourile în afara cercului unitar şi sunt

prime între ele.

Există o mare varietate de variabile instrumentale

pentru cazuri particulare de forme ale filtrului K(q-1) şi

polinoamelor C(q-1) şi D(q-1).

Pentru fiecare din aceste cazuri se determină parametrii astfel

ca relaţiile (5.78), (5.79) să fie îndeplinite. Dacă de exemplu

nc=nd=na relaţia (5.80) devine:

16.05.2016

3

TT nbnakukukuqD

qKDCSnbkuqD

kuqDnakuqCkuqCqD

qKkz

)(),...,2(),1()(

)(),()]()(...

)...1()(),()(),...,1()([)(

)()(

1

11

111

1

1

(5.82)

unde S(-C, D) este matricea Sylvester asociată polinoamelor -C şi D :

ba

b

nn

a

na

nana

na

b

na

nana

n

na

n

dd

dd

d

n

ccc

ccc

ccc

DCS

... 1 ... 0 0

..........................................................................

0 0 .. ...d 1 0

0 0 0 ...d 1

... 0 0

0 0 0.... ... 0

0 0 .... 0 0 ...

),(

1

11

1

10

10

10

Dacă C(q-1) şi D(q-1) sunt polinoame prime între ele rezultă

că rangul matricei S(-C, D) este egal cu suma na + nb,

deci matricea S este nesingulară.

În cazul în care polinoamele C şi D au k zerouri comune, se arată în [28] că rang S(-C, D)= na + nb -k. In relaţia (5.82) matricea S(-C, D) reprezintă o transformare

liniară a vectorului de variabilă instrumentală, ceea ce nu

afectează estimatorul. Ca urmare, vectorul z(k) poate fi

TnbnakukuqD

qKkz )(..).........1(

)(

)()(

1

1

(5.83)

ceea ce înseamnă că de fapt polinomul C(q-1) nu afectează

estimatorul.

În particular, dacă K(q-1) = D(q-1) se găseşte varianta

Tba nnkukukz )(..).........1()( (5.84)

O altă variantă de alegere a vectorului de variabilă instrumentală,

în care valorile ieşirii sunt întârziate cu i intervale de

eşantionare, ceea ce slăbeşte corelaţia cu perturbaţia, este

definită prin

16.05.2016

4

Tba nkukunikyikykz )(),....,1(),(),....,1()( (5.85)

Condiţia a doua de consistenţă este îndeplinită dacă perturbaţia

v(k) este de medie alunecătoare de ordin . i

.Dacă sistemul este descris de ecuaţia

)()(,)( 1 kvkuqHky (5.86)

atunci )(,)( 1 kuqHkx

ieşirii ( neafectată de zgomot) a sistemului. Metoda de variabilă

instrumentală idealizată se bazează pe alegerea

este componenta deterministă a

Tba nkukunkxkxkkz )(),...,1(),(),....,1()(~)( (5.87)

În ipoteza unui semnal de intrare persistent o astfel de alegere

asigură îndeplinirea ambelor condiţii de consistenţă. Într-adevăr,

dacă se notează

Tnakvkvkv 0,0,...0),(),....,1()(~ (5.87)

rezultă imediat că Tnbkukunakykyk )](),...,1(),(),...,1([)(

)(~)(~)( kvkk

)()(~0

0)(~)(~)()(

kvkE

kkEkkzER TT

şi (5.89)

În ipoteza unui semnal de intrare persistent o astfel de alegere

asigură îndeplinirea ambelor condiţii de consistenţă. Într-

adevăr, dacă se notează

Tnakvkvkv 0,0,...0),(),....,1()(~ (5.88)

deoarece se presupune că intrarea şi perturbaţia sunt necorelate.

Matricea R este simetrică şi cel puţin nenegativ definită şi

este pozitiv definită dacă intrarea este semnal persistent.

Această variantă prezintă interes teoretic, nu şi practic, deoarece este necunoscut. Cunoaşterea lui )(~ k )(~ k

ar însemna cunoaşterea modelului părţii deterministe

a sistemului, întrucât x(k) este intrarea

16.05.2016

5

filtrată de partea deterministă.

Cu toate acestea, se poate imagina un algoritm iterativ în care

modelul adevărat, reprezentat prin parametrii θ*, necesar

pentru determinarea variabilelor instrumentale idealizate, este

înlocuit cu o aproximaţie a acestuia.

),(~)ˆ,( kkz

)()(ˆ 1 ikuqH )(ˆ 1qH )(1

qH

Fie în care ieşirile libere de zgomot x(k-i) sunt

, cu un estimator al lui

, corespunzând vectorului parametrilor .

.

înlocuite prin

Algoritmul iterativ se bazează pe următoarea relaţie de recurenţă.

N

k

kN

k

Tkk kykzN

kkzN 1

1

1

1 )(ˆ,1

)(ˆ,1ˆ (5.90)

cu iniţializarea 0 , care poate fi de exemplu, un estimator CMMP

5.5.3. Estimarea parametrilor prin metoda variabilelor

instrumentale iv4 Pentru estimarea parametrilor prin metoda variabilelor

instrumentale în MATLAB [85], [87] se aplică

funcţia iv4.

Algoritmul cuprinde mai multe etape: în prima etapă se aplică

funcţia arx şi se obţine un model folosit la generarea variabilelor

instrumentale; în etapa a doua se aplică metoda variabilelor

instrumentale iv4; reziduurile obţinute în etapa a doua sunt

modelate printr-un model AR în etapa a treia; în etapa a patra se

aplică din nou metoda variabilelor instrumentale iv4

Exemplul 5.4. Se consideră sistemul din exemplul 5.2, cu un

răspuns similar celui din fig. 5.3. Să se determine un model cu

structura nn = [2 2 1] prin metoda iv4, utilizând programul

pismnr.m din Anexa 3. Ştiind structura nn şi datele intrare-

ieşire din matricea z, cu funcţia iv4 se obţine modelul thi, care

este afişat cu funcţia present

>>thi=iv4(z,nn)

>>present(thi)

Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t)

16.05.2016

6

A(q) = 1 - 1.199 (+-0.00315) q^-1 + 0.5999 (+-0.002534) q^-2

B(q) = 0.9942 (+-0.005995) q^-1 + 0.4983 (+-0.00876) q^-2

Estimated using IV4 from data set z

Loss function 0.0107966 and FPE 0.0110884 Se determină reziduurile ei ale modelului IV4 cu funcţia resid şi

se reprezintă grafic în fig.5.9.

>>ei=resid(z,thi);

0 5 10 15 20 25-1

-0.5

0

0.5

1Functia de corelatie a rezidualilor. Iesirea y- metoda IV

-30 -20 -10 0 10 20 30-0.2

-0.1

0

0.1

0.2Functia de intercorelatie intre intrarea u si rezidualii din iesirea y

lag

Fig. 5.9

Răspunsul sistemului, y şi a modelului IV4, ym, pentru primele 50

de eşantioane, se reprezintă grafic în fig.5.10.

>>plot(t1,y(1:50),'b',t1,ym(1:50),'g');grid

>>title('Răspuns sistem y şi răspuns model – ym. Metoda IV4')

0 10 20 30 40 50-6

-4

-2

0

2

4

6Raspuns sistem y si raspuns model IV4 ym

y

ym

Fig. 5.10 5.6. Metode recursive de identificare parametrică a

sistemelor Metodele de identificare a parametrilor de tip off-line

prezentate anterior prelucrează simultan întregul şir de

măsurători intrare-ieşire.

16.05.2016

7

În cazul metodelor de identificare recursive sau on-line

datele sunt prelucrate pe măsură ce ele devin disponibile prin

măsurători.

Procedurile on-line se impun în următoarele situaţii

[28],[38],[76]:

- când se doreşte continuarea experimentului de identificare

până la obţinerea unei precizii aprioric impuse;

- când sistemul are parametri variabili în timp;

- în cazul sistemelor de control adaptiv conform schemei bloc

din fig. 5.11.

Parametrii procesului fiind variabili în timp, procesul nu

poate fi reprezentat printr-un model valabil pentru orice

moment de timp. Se utilizează un model variabil în timp pentru a determina

parametrii unui regulator, de asemenea, variabili în timp. În

acest fel regulatorul va fi dependent de comportarea previzibilă

a procesului.

)1(ˆ N

)(ˆ N

Algoritmii de estimare on-line, pot fi relativ simplu obţinuţi din cei

off-line printr-o organizare corespunzătoare a calculelor, astfel ca

estimaţia bazată pe N+1 date, să poată fi calculată

şi alte câteva măsurători anterioare. Estimaţia recursivă este

totuşi o aproximare a estimaţiei off-line şi deci este mai puţin

precisă.

recursiv utilizând numai , măsurătorile la momentul N+1

Identificare

recursivă

Proces (sistem)

Regulator

u

u

v(t)

perturbaţie

Parametri estimaţi

y (ieşire)

Fig. 5.11

Între algoritmii de estimare recursivi (on-line) şi versiunile lor

în timp real trebuie făcută o distincţie. O procedură

de estimare va fi de tip real dacă este o

16.05.2016

8

procedură recursivă şi în plus, estimatorul parametrilor poate

urmări variaţiile în timp ale parametrilor sistemului.

Metodele de identificare on-line cele mai utilizate sunt:

1. Metoda recursivă a celor mai mici pătrate (RCMMP).

2. Metoda recursivă a celor mai mici pătrate generalizate

(RCMMPG)

3. Metoda recursivă a variabilelor instrumentale (RVI)

4. Metoda recursivă a verosimilităţii maxime (RVM) 5.6.1. Metoda recursivă a celor mai mici pătrate (RCMMP)

5.6.1.1. Principiul metodei recursive a celor mai mici pătrate

Pentru un sistem dinamic liniar monovariabil discret descris de

ecuaţia stocastică cu diferenţe

)()( ; )()()()()( 11 kekvkvkuqBkyqA

***

1*1

1

***

1*1

1

......)(

......1)(

nbnb

nana

qbqbqB

qaqaqA

(5.91)

se consideră un model de forma

nbna

nana

qbqbqbqB

qaqaqA

kuqBkyqA

.....)(

......1)(

)()()()(

22

11

1

11

1

11

(5.92)

Se definesc vectorii : (vectorul parametrilor modelului) şi (k) :

Tnbna bbaaa ],...,,...,,[ 121

Tnbkukunakykyk )]()...1(),(),...,1([)( (5.93)

Modelul (5.92) se scrie sub forma

)()( kky Tm (5.94)

Pentru o funcţie criteriu

N

k

TN

k

N

km

kkyN

kuqBkyqAN

kykyqAN

V

1

2

1

211

1

212

])()([1

)]()()()([1

)]()()[(1

)(

(5.95)

16.05.2016

9

se determină un estimator al celor mai mici pătrate pentru N

măsurători

N

k

N

k

T

N

k

N

k

T

kykkk

kykN

kkN

N

1

1

1

1

1

1

)()())()((

)()(1

))()(1

()(ˆ

(5.96)

Pentru o funcţie criteriu generală

N

k

TkN

N

km

kN

kkyN

kykyqAN

V

1

2

1

212

])()([1

)]()()[(1

)(

(5.97)

unde λ este un coeficient de ponderare a erorii pătratice,

pentru N măsurători, se obţine estimatorul

N

k

kNN

k

TkN kykkkN1

1

1

)()())()(()(ˆ (5.98)

Se introduce notaţia 1

1

])()([)(

N

k

TkN kkNP (5.99)

şi estimatorul (5.98) devine

N

k

kN kykNPN1

)()()()(ˆ (5.100)

Relaţia (5.99) poate fi scrisă într-o formă recursivă

11

1

1

11

1

1

)]1()1()([

)]1()1()()([

])()([)1(

NNNP

NNkk

kkNP

T

TN

k

TkN

N

k

TkN

(5.101)

Se ţine seama de următoarea identitate de inversiune

matricială [76]

111111 )1()( AbbAbbAAbbA TTT (5.102)

)(1 NPA )1( Nb

A şi b fiind matrici de dimensiuni corespunzătoare. Această

identitate se verifică uşor prin calcul direct.

Utilizând relaţia (5.102), pentru şi

, din relaţia (5.101) se obţine

16.05.2016

10

])1()()1(

)()1()1()()([

1)()1(

)]1()(

)1(1)[1()()(

)1( 1

NNPN

NPNNNPNP

NPN

NNP

NNNPNP

NP

T

TT

T

(5.103)

În relaţia (5.103) s-a avut în vedere că )1()()1( NNPNT

este o mărime scalară. Cu notaţia

)1()()1(

)1()()1(

NNPN

NNPNK

T

(5.104)

relaţia (5.103) se scrie în forma

)()]1()1([1

)]()1()1()([1

)1(

NPNNKI

NPNNKNPNP

T

T

(5.105)

Pe baza relaţiilor (5.100)(5.105) se poate determina

valoarea estimatorului )1(ˆ N

:

)1()1()(ˆ

](ˆ)1()1()[1()(ˆ

)1()1()(ˆ)1()1()(ˆ

)1()]1()()1()1()1()([1

)(ˆ)1()1()(ˆ

)1()1()()1()1(

)1()1()()()(

)()1()1()()()(

])1()1()()()[1(

)()()1()1(ˆ

1

1

1

1

1

1

NNKN

NNNyNKN

NyNKNNNKN

NyNNPNNKNNP

NNNKN

NyNNPNNK

NyNNPkyk

NPNNKkykNP

NyNkykNP

kykNPN

T

T

T

T

T

N

k

kN

N

k

TkN

N

k

kN

N

k

kN

(5.106)

S-a ţinut seama că din (5.104) rezultă

16.05.2016

11

)1()1()()1()1()1()( NKNNPNNKNNP T

De asemenea în (5.106) s-a introdus notaţia

)(ˆ)1()1()1( NNNyN T (5.107)

Relaţiile (5.100) (5.107) permit calculul recursiv (on-line) al).

estimatorului celor mai mici pătrate (CMMP Teorema 5.4. Estimatorul celor mai mici pătrate satisface

următoarele ecuaţii recursive:

)1()1()(ˆ)1(ˆ kkKkk

)(ˆ)1()1()1( kkkyk T

)1()()1(

)1()()1(

kkPk

kkPkK

T

)]()1()1()([1

)1( kPkkKkPkP T

(5.108)

(5.109)

(5.110)

(5.111)

Demonstraţie. Relaţiile (5.109) (5.111) se obţin din relaţiile.

(5.100)- (5.107) dacă se notează : x(N)=x(k) ,

x(N+1)=x(k+1), unde x(.) sunt mărimile care intervin în aceste relaţii

Se consideră că v(k)=e(k) este zgomot alb şi că u(k) este un

semnal persistent de ordin bn

.

)(ˆ k

Nk

Deoarece prin construcţie este estimatorul CMMP (off-line)

date, rezultă că bazat pe

c.p.1 )(ˆlim

k

k(5.112)

θ* este vectorul parametrilor adevăraţi ai sistemului.

)1(ˆ k

)(ˆ k

)1( k

Egalitatea (5.112) rezultă din consistenţa estimatorului CMMP

off-line. Convergenţa estimatorului RCMMP este deci simplu

de stabilit.

Relaţia (5.108) are un caracter intuitiv, noua estimaţie

se obţine aplicând estimaţiei vechi o corecţie proporţională

între valoarea măsurată a ieşirii y(k+1)

.

cu diferenţa

la momentul k+1 şi valoarea predictată a acesteia,

16.05.2016

12

)(ˆ)1( kkT , bazată pe modelul reprezintă eroarea de predicţie de pas optimală.

)(ˆ k

)1( k

)1( k

Coeficientul de corecţie K(k+1) este modificat la rândul său la

fiecare pas.Relaţiile (5.108) - (5.111) pot fi reduse la două dacă

. Din relaţia (5.110) rezultă: se elimină K(k+1) şi

)]1()()1()1()1()([1

)1( kkPkkKkkPkK T

(5.113)

Înmulţind cu φ(k+1) relaţia (5.111) se obţine

)1()1(1()1()1(1

)]1()()1()1()1()([1

)1()1(

kkPkKkKkK

kkPkkKkkPkkP T

(5.114)

Înlocuind K(k+1) din (5.112) în (5.109) şi explicitând P(k+1)

algoritmul recursiv se reduce numai la două relaţii :

)]()1()1()1()([1

)1(

)](ˆ)1()1()[1()1()(ˆ)1(ˆ

kPkkkPkPkP

kkkykkPkk

T

T

(5.115)

1)()1()1()(

)1(

kPkkI

kPkP

T

(5.115)

Aceasta justifică iniţializarea algoritmului recursiv cu valorile

)0(

Trebuie menţionat că matricea K(k) are semnificaţie fizică, fiind

matricea Kalman corespunzătoare modelului exprimat în ecuaţii

de stare, pentru care s-a dedus estimatorul celor mai mici pătrate.

Pentru N eşantioane estimatorul CMMP se poate scrie şi în

forma:

şi P(0).

T

T

T

TT

bnkukuankykyk

NyyyNy

NN

NyNNNN

)]ˆ(),.....1(),ˆ(),....,1([)(

)](.,),........2(),1([)(

)](..,),........2(),1([(

)()()]()([)(ˆ 1

(5.116)

Sistemul este descris de ecuaţia

)()()( * NEeNNy (5.117)

16.05.2016

13

unde Te NeeNE )](,),........1([)(

cu matricea de covarianţă λ*2I. Înlocuind y(N) din (5.117) în

(5.116) rezultă

este secvenţa de zgomot alb

1

*

1*

)]()([)(

)()()(

)()()]()([)(ˆ

NNNP

NENNP

NENNNN

T

eT

eTT

(5.118)

Matricea de covarianţă a estimatorului este:

)()()()()(

)()()()()()(

)]()()()()()([

]))(ˆ)()(ˆ[()(ˆcov

2*2*

**

NPNPINNNP

NPNNENEENNP

NPNNENENNPE

NNEN

T

TTee

T

TTee

T

T

(5.119)

Din (5.117) rezultă că P(N) este proporţională cu matricea de

covarianţă a estimatorului . Pentru că dispersia λ*2 nu este )(ˆ N

în general cunoscută, ea poate fi înlocuită cu o estimaţie

bazată pe N date.

N

k

T NVNkkyN

N1

22 )())(ˆ)()((1

)(ˆ (5.120)

Se poate demonstra că V(N) satisface relaţia de recurenţă

următoare

)1()()1(

)1()()1(

2

NNPN

NNVNV

T

(5.121)

5.6.1.2. Consideraţii numerice

În algoritmul RCMMP efectul acumulării erorilor de rotunjire

asupra algoritmului poate să fie destul de important. Pentru

reducerea acestui efect se ţine seama că relaţia (5.111) se

poate scrie şi sub forma

)1()1(1

)]1()1()[()]1()1([1

)1(

kKkK

kkKIkPkkKIkP

T

TTT

(5.122)

Noua expresie (5.122) pentru P(k+1) are următoarele

avantaje faţă de expresia (5.111) :