Subspace Identification Methods

T´ ema 7

Subspace Identification Methods

OaF - t´ ema 7


Co jsou to subspace metody? Metody pro identifikaci parametr˚ u lineárn´ıch stavov´ ych model˚ u z experimentáln´ıch vstupnˇe/v´ ystupn´ıch dat. • zaloˇzeny na geometrick´ ych operac´ıch (projekce, pr˚ uniky) a numerické lineárn´ıch algebˇre • alternativa k rozˇs´ıˇren´ ym prediction error metodám PEM (napˇr. ML identifikace ARX model˚ u pomoc´ı nejmenˇs´ıch ˇctverc˚ u) Vlastnosti: • stejn´ a sloˇzitost identifikace pro systémy s jedn´ım vstupem a v´ ystupem (SISO) a pro systémy s v´ıce vstupy a v´ ystupy (MIMO) • mal´ y poˇcet uˇzivatelem nastavovan´ ych parametr˚ u (prakticky pouze ˇrád systému, jehoˇz odhad dokonce i samy poskytuj´ı) • numericky robustn´ı implementace pomoc´ı SVD a QR faktorizac´ı • implicitn´ı redukce ˇr´ adu modelu • nevhodné pro malé soubory dat Oznaˇcované jako 4SID metody (Subspace State-Space System IDentification). N´ azev ”subspace” vyjadˇruje skuteˇcnost, ˇze parametry stavového modelu mohou b´ yt z´ıskány z ˇr´ adkov´ ych a sloupcov´ ych prostor˚ u urˇcit´ ych matic vytvoˇren´ ych ze vstupnˇe/v´ ystupn´ıch dat.

c P. Trnka, 2007 °

2

OaF - t´ ema 7

7.1


Motivace

Identifikace syst´ emu v´ıce vstupy a v´ıce v´ ystupy (MIMO) • ARX, ARMAX, OE modely vyˇzaduj´ı správnou strukturáln´ı parametrizaci, tj. správnou volbu ˇrádu jednotliv´ ych polynom˚ u. Napˇr´ıklad pro MIMO ARX model s m vstupy a ` v´ ystupy A(q)y(t) = B(q)u(t) + e(t),  a (q) · · · a1m (q)  11 .. ..  .. A(q) =  . . .  a`1 (q) · · · a`m (q)

   , 

   B(q) =  

 b11 (q) · · · .. .. . .

b1m (q) .. .

b`1 (q)

b`m (q)

···

  , 

aij (q) a bij (q) jsou polynomy jejichˇz jednotlivé ˇrády je nutno zvolit. Tato volba nen´ı zˇrejmá, avˇsak v´ yznamnˇe ovlivˇ nuje v´ ysledn´ y identifikovan´ y model. • Z hlediska parametrizace MIMO modelu je snaˇzˇs´ı pouˇz´ıt stavov´ y model, kde jedin´ ym strukturn´ım parametrem je ˇ r´ ad syst´ emu n. Predikˇ cn´ı schopnosti modelu • PEM metody hledaj´ı modely s jednokrokovˇe optimáln´ımi predikcemi • 4SID metody hledaj´ı modely optimalizované na v´ıce krokové predikce - vhodné pro aplikace jako MPC Pozn.: ARX nen´ı identifikaˇcn´ı metoda, ale struktura modelu. Metoda odhadu parametr˚ u pomoc´ı nejmenˇs´ıch ˇctverc˚ u (LS) pak vych´ az´ı z principu maximáln´ı vˇerohodnosti (ML). c P. Trnka, 2007 °

3

OaF - t´ ema 7


Uhel mezi vektory signalu u(t) a e(t) v zavislosti na poctu vzorku 180 Stredni hodnota Standardni odchylka

160 140 120 Uhel [deg]

Obecn´ y princip 4SID metod Podstatou subspace metod zaloˇzen´ ych na geometrick´ ych vlastnostech prostor˚ u generovan´ ych signály vstup˚ u, v´ ystup˚ u a ˇsum˚ u je ortogonalita mezi nez´ avisl´ ymi sign´ aly. Máme-li dva nezávislé skal´ arn´ı signály u(t) a e(t) (E {e(t)} = 0), napˇr´ıklad diskrétn´ı signály   v(t) u(t) = 10 + 10 sin(t/2) + v(t),  ∼ N (0, I).  w(t) e(t) = w(t),

100 80 60 40 20 0

a posklád´ ame-li N vzork˚ u z tˇehto signál˚ u do vektor˚ uua v, potom u ´hel α mezi tˇemito vektory cos(α) =

−20 0 10

1

10

2

10 Pocet vzorku

3

10

4

10

uT v |u||v|

se bude bl´ıˇzit π/2. Graf ukazuje stˇredn´ı hodnotu a smˇerodatnou odchylku u ´hlu ze 100 realizac´ı pro kaˇzdou délku N . Demonstruje, ˇze i pro malé délky dat je ve stˇredn´ı hodnotˇe u ´hel π/2, avˇsak s velkou odchylkou pro jednotlivé realizace. S rostouc´ım poˇctem vzork˚ u se odchylka rychle zmenˇsuje.

c P. Trnka, 2007 °

4

OaF - t´ ema 7

7.2


Geometrick´ e a algebraick´ e n´ astroje

adkové maticové prostory • sloupcové / ˇr´ • Hankelovy matice • kolm´ a projekce (orthogonal projection) • ˇsikm´ a projekce (oblique projection) aln´ı u ´hly a smˇery (principal angles and directions) • principi´ • singul´ arn´ı rozklad (Singular Value Decomposition - SVD) ˇ adkov´ R´ y prostor matice ˇ adkov´ R´ y prostor matice A ∈ Rm×n , oznaˇcovan´ y jako row(A), je podprostor Rn tvoˇren´ y vˇsemi lineárn´ımi kombinacemi (lineárn´ım obalem) ˇr´ adkov´ ych vektor˚ u matice A row(A) = ha1 , . . . , am i ,

kde A = [ aT1

T aTm ] ,

...

symbol h•i znaˇc´ı lineárn´ı obal v nˇem uzavˇren´ ych vektor˚ u •. Sloupcov´ y prostor matice Sloupcov´ y prostor matice A ∈ Rm×n , oznaˇcovan´ y jako col(A), je podprostor Rm tvoˇren´ y vˇsemi lineárn´ımi kombinacemi (lineárn´ım obalem) sloupcov´ ych vektor˚ u matice A col(A) = ha1 , . . . , an i ,

kde A = [ a1

...

an ].

c P. Trnka, 2007 °

5

OaF - t´ ema 7


Hankelovy matice Matice se stejn´ ymi prvky na vedlejˇs´ıch diagonálách. Z posloupnosti {a(t)}m+n−1 lze Hankelovu matici 1 m×n H∈R vytvoˇrit jako  a(1) a(2) a(3) . . .   a(2) a(3)   ..  . H =  a(3)   .. ..  . .  a(m) ...

 a(n) .. .

Hankelova matice

         

1 2 3 4 5 6 7

a(m + n − 1)

2

4

6

8

10

12

14

16

18

20

kde hodnota prvku na pozici (i, j) závis´ı pouze na souˇctu i + j − 1 hi,j = a(i + j − 1). Jednotlivé prvky mohou b´ yt také vektory nebo matice, potom se jedná o blokovou Hankelovu matici. Proˇ c Hankelovy matice? Je-li {y(t)}k0 odezva systému n-tého ˇr´ adu s nulov´ ymi poˇcáteˇcn´ımi podm´ınkami x(0) na libovoln´ y vstupn´ı signál, potom Hankelova matice Y vznikl´ a z této odezvy má hodnost nanejv´ yˇse n. Nav´ıc je-li systém dostateˇcnˇe vybuzen, potom sloupcov´ y prostor col(Y ) je invariantn´ı a je shodn´ y se sloupcov´ ym prostorem rozˇs´ıˇrené matice pozorovatelnosti (viz. rovnice rozˇs´ıˇreného stavového modelu). c P. Trnka, 2007 °

6

OaF - t´ ema 7


SVD - Singular Value Decomposition Singulárn´ı rozklad faktorizuje obdéln´ıkovou matici A ∈ Rm×n na souˇcin tˇr´ı matic     h i S1 0 h iT σ1 .  V1 V2 , .. A = U S V T = U1 U2  , S1 =  0 0 σ` kde S ∈ Rm×n je diagon´ aln´ı se singulárn´ımi ˇc´ısly na diagonále a U ∈ Rm×m resp. V ∈ Rn×n jsou ortonorm´ aln´ı se sloupci tvoˇren´ ymi lev´ ymi resp. prav´ ymi singulárn´ımi vektory, S1 ∈ R`×` obsahuje nenulová singulárn´ı ˇc´ısla v nerostouc´ım poˇrad´ı σ1 ≥ σ2 ≥ · · · ≥ σ` . Siln´ y nástroj lineárn´ı algebry s ˇsirok´ ym vyuˇzit´ım: • Urˇcen´ı hodnosti rank(A) = `, posouzen´ı ”bl´ızkosti k singularitˇe” z nejmenˇs´ıho singulárn´ıho ˇc´ısla. • Maticov´ a pseudoinverze A+ = V S −1 U T (S −1 má na diagonále inverze nenulov´ ych singulárn´ıch ˇc´ısel). ˇ sen´ı u • Reˇ ´lohy nejmenˇs´ıch ˇctverc˚ u Zθ = b i v pˇr´ıpadˇe regresoru Z bez plné sloupcové hodnosti jako θˆ = Z + b a posouzen´ı numerické podm´ınˇenosti ˇreˇsen´ı (citlivosti na zmˇeny v datech) podle pod´ılu nejvˇetˇs´ıho a nejmenˇs´ıho singulárn´ıho ˇc´ısla (condition number κ = σmax /σmin ). • Nalezen´ı ortonormáln´ı báze prostoru j´ adra A tvoˇrené prav´ ymi singulárnimi vektory V2 odpov´ıdaj´ıc´ı nulov´ ym singulárn´ım ˇc´ısl˚ um. • Nalezen´ı ortonormáln´ı báze sloupcov´ eho prostoru A tvoˇrené lev´ ymi singulárnimi vektory U1 odpov´ıdaj´ıc´ı nenulov´ ym singulárn´ım ˇc´ısl˚ um. • Aproximace matic´ı niˇ zˇ s´ı hodnosti - nalezen´ı matice A0 s niˇzˇs´ı hodnost´ı k, která je nejlepˇs´ı aproximac´ı A ve smyslu minimalizace Frobeniovy maticové normy 2

min kA − A0 kF A0

za podm´ınky rank(A0 ) = k,

c P. Trnka, 2007 °

7

OaF - t´ ema 7


Frobeniova norma k•kF je odmocnina ze sumy kvadrát˚ u jednotliv´ ych prvk˚ u matice. Aproximaci A0 z´ıskáme zanedb´ an´ım singulárn´ıch vektor˚ u odpov´ıdaj´ıc´ı mal´ ym singulárn´ım ˇc´ısl˚ um. Vyuˇzit´ım Matlabovské notace A0 = U (:, 1 : k)S(1 : k, 1 : k)V (:, 1 : k)T . Pˇ r´ıklad: Ztr´ atov´ a komprese obrázk˚ u – matice A obsahuje obrázek o rozmˇerech 512 × 512 = 256 kB – SVD faktorizac´ı z´ısk´ ame tˇri matice U , S, V T se stejn´ ymi rozmˇery 512 × 512, avˇsak pro ztrátovou kompresi obrázku staˇc´ı ukládat k lev´ ych a prav´ ych singulárn´ıch vektor˚ u pˇr´ısluˇsn´ ych k nejvˇetˇs´ım singulárn´ım ˇc´ısl˚ um – napˇr´ıklad pro k = 50 uloˇz´ıme D 1 = U (:, 1 : k)S(1 : k, 1 : k)1/2 a D 2 = S(1 : k, 1 : k)1/2 V (:, 1 : k)T o celkové velikosti 50 kB – zpˇetn´ a rekonstrukce A0 = D 1 D 2 (matice hodnosti k) k=10 (96.1% compression)

c P. Trnka, 2007 °

k=50 (80.5% compression)

Original

8

OaF - t´ ema 7

7.3


Stavov´ a realizace

−1 Metoda autor˚ u Ho & Kalman (1966) umoˇzn ˇuje pro danou impulsn´ı posloupnost {h(t)}N naj´ıt parametry 0 A, B, C, D odpov´ıdaj´ıc´ıho stavového modelu. Bez pˇr´ıtomnosti ˇsumu, pro systém n-tého ˇrádu, staˇc´ı znát alespoˇ n 2n prvk˚ u posloupnosti.

Zaloˇzena na moˇznosti napsat Hankelovu matici impulsn´ıch posloupnost´ı jako souˇcin rozˇs´ıˇrené matice pozorovatelnosti Γi a rozˇs´ıˇrené matice ˇriditelnosti ∆i   C     h i  CA  , ∆i = B AB . . . Ai−1 B , Γi =  ..     .   CAi−1 kde i je vˇetˇs´ı neˇz ˇr´ ad systému n. Impulsn´ı odezva stavového modelu je dána parametry A, B, C, D jako h(0) = Cx(0) + Du(0) = D h(1) = Cx(1) + Du(1) = C(Ax(0) + Bu(0)) + Du(1) = CB h(2)

= Cx(2) + Du(2) = C(Ax(1) + Bu(1)) + Du(1) = C(A(Ax(0) + Bu(0)) + Bu(1)) + Du(1) = CAB .. .

h(i) = CAi−1 B

c P. Trnka, 2007 °

9

OaF - t´ ema 7


Uspoˇra´d´ an´ı impulsn´ı posloupnost do Hankelovy matice    h(1) h(2) CB CAB       . .. .. H =  h(2)  =  CAB .    h(N − 1)

   ,  CAN −2 B

lze napsat jako souˇcin H = Γi ∆i . Pro pozorovateln´ y a ˇriditeln´ y systém ˇr´ adu n maj´ı rozˇs´ıˇrené matice pozorovatelnosti a ˇriditelnosti plnou sloupcovou resp. ˇ r´ adkovou hodnost ⇒ rank(H) = n. Toho lze vyuˇz´ıt a singulárn´ım rozkladem H = U SV T , tak dostaneme n nenulov´ ych singulárn´ıch ˇc´ısel, jejichˇz poˇcet urˇ cuje ˇ r´ ad syst´ emu. Rozˇs´ıˇrenou matici ˇriditelnosti a pozorovatelnosti lze spoˇc´ıtat jako Γi = U (:, 1 : n)S(1 : n, 1 : n)1/2 ,

∆i = S(1 : n, 1 : n)1/2 V (:, 1 : n)T ,

kde odmocnina z S(1 : n, 1 : n) je snadná, nebot’ se jedná o diagonáln´ı matici. Nalezené Γi a ∆i nemus´ı ˇc´ıselnˇe odpov´ıdat p˚ uvodn´ımu systému, coˇz je dáno neurˇcitou volbou báze stavového modelu. Urˇ cen´ı parametr˚ u A, B, C, D ze znalosti Γi a ∆i Pro odhad parametr˚ u stavového modelu lze vyuˇz´ıt napˇr´ıklad tzv. invariantnost Γi v˚ uˇci posunu Γi = Γi A, c P. Trnka, 2007 °

10

OaF - t´ ema 7


kde Γi resp. Γi jsou rozˇs´ıˇrené matice pozorovatelnosti Γi bez prvn´ıho resp. bez posledn´ıho blokového ˇrádku     CA C        2   CA   CA  ,  . Γ Γi =  = i .. ..         . .     i−1 i−2 CA CA Matici A tak lze z´ıskat napˇr´ıklad pomoc´ı nejmenˇs´ıch ˇctverc˚ u, matici B jako prvn´ı blokov´ y sloupec ∆i , matici C jako prvn´ı blokov´ y ˇr´ adek Γi a matici D jako prvn´ı ˇclen impulsn´ı posloupnosti +

A

=

Γi Γi ,

B

=

∆i (:, 1 : m),

C

=

Γi (1 : l, :),

D

=

h(0),

kde m je poˇcet vstup˚ u a l je poˇcet v´ ystup˚ u. Realizace z ˇ sumem zat´ıˇ zen´ e impulsn´ı odezvy • Nenulov´ ych singulárn´ıch ˇc´ısel bude v´ıce jak n. • Pro n´ızkou m´ıru ˇsumu bude v posloupnosti singulárn´ıch ˇc´ısel patrné skokové zmenˇsen´ı pouˇzitelné pro odhad ˇr´ adu systému. Patrnost skoku se bude s rostouc´ı m´ırou ˇsumu rychle vytrácet. • Pro zvolen´ y ˇr´ ad poskytne algoritmus takov´ y stavov´ y model, jehoˇz impulsn´ı posloupnost bude m´ıt pˇribliˇznˇe minim´ aln´ı souˇcet kvadr´ at˚ u chyb odchylek od zaˇsumˇené impulsn´ı posloupnosti. c P. Trnka, 2007 °

11

OaF - t´ ema 7


−3

4

Pˇ r´ıklad - pokles singulárn´ıch ˇc´ısel pˇri realizaci ze zaˇsumˇené impulsn´ı odezvy • systém P (s) =

1 (s+1)(s+3)(s+5)

Impulsni odezvy

x 10

Skutecna odezva Zasumena SNR~20dB Zasumena SNR~40dB

3.5 3 2.5

vzorkovan´ y s periodou Ts = 0.1s

2

39

• k dispozici 40 vzork˚ u impulsn´ı posloupnosti {h(t)}0

1.5 1

• zat´ıˇzen´ı ˇsumem s pomˇerem signál k ˇsumu 0dB, 20dB a 40dB

0.5 0

0

5

10

15

20

25

30

35

40

45

Singularni cisla pro ruznou miru sumu SNR=0dB

0

SNR~20dB

0

10

SNR~40dB

0

10

10

−5

10

−2

−2

10

10

−10

10

−4

−4

10

−15

10

10

−20

10

−6

0

5

10

15

20

10

−6

0

5

10

15

20

10

0

5

10

15

20

Shrnut´ı • Pro systémy s v´ıce vstupy a v´ ystupy funguje realizaˇcn´ı metoda zcela stejnˇ e (h(i) jsou matice). • Nutnost znát nebo z´ıskat impulsn´ı posloupnost je znaˇcnˇe omezuj´ıc´ı. Tuto nev´ yhodu odstraˇ nuj´ı právˇe 4SID metody, které pracuj´ı obdobnˇe, avˇsak vystaˇc´ı si se vstupnˇe/v´ ystupn´ımi experimentáln´ımi daty. c P. Trnka, 2007 °

12

OaF - t´ ema 7

7.4


Rozˇ s´ıˇ ren´ y stavov´ y model

Stavov´ y model v inovaˇ cn´ım tvaru Z hlediska identifikace je v´ yhodné pouˇz´ıvat stochastick´ y stavov´ y model v inovaˇ cn´ım tvaru x(t + 1)

=

y(t) =

Ax(t) + Bu(t) + Ke(t),

cov {e(t)} = Re ,

Cx(t) + Du(t) + e(t),

kter´ y z´ısk´ ame aplikac´ı Kalmanova filtru na stochastick´ y stavov´ y model v obvyklém tvaru s ˇsumem procesu v(t) a ˇsumem mˇeˇren´ı e0 (t)      v(t)  x(t + 1) = Ax(t) + Bu(t) + v(t), Q 0  = . cov   e0 (t)  y(t) = Cx(t) + Du(t) + e0 (t), 0 R Kalmanovo zes´ılen´ı K z´ısk´ ame z algebraické Riccatiho rovnice ³ ´−1 T K = AP C CP C T + R , ³ ´ P = AP AT − K C T P C + R K T + Q, a kovarianci Re Re = CP C T + R. • oba modely jsou ekvivalentn´ı z hlediska stˇredn´ı hodnoty a kovariance v´ ystupu y(t) • inovaˇcn´ı model má oproti obvyklému modelu menˇs´ı poˇcet stupˇ n˚ u volnosti ve stochastickém subsystému, které mohou byt jednoznaˇcnˇeji identifikovány ze vstupnˇe/v´ ystupn´ıch dat. c P. Trnka, 2007 °

13

OaF - t´ ema 7


Rozˇ s´ıˇ ren´ y stavov´ y model Známe-li stav x(t), posloupnost i vstup˚ u {u(t), u(t + 1), . . . , u(t + i − 1)} a inovac´ı {e(t), e(t + 1), . . . , e(t + i − 1)}, potom pro v´ ystup y(t) na 0 aˇz i − 1 krok˚ u m˚ uˇzeme psát         y(t) C D u(t) I                 D I  y(t+1)   CA   CB  u(t+1)   CK  =    + .. ..    ..  x(t)+  ..    .. ..    .   .    . . . CB . CK         y(t+i−1)

CA

i−1

i−2

CA B

...

CB

D u(t+i−1)

i−2

CA K

Pro posloupnost poˇc´ ateˇcn´ıch stav˚ u {x(t), x(t + 1), . . . , x(t + j − 1)}   y(t) y(t+1) . . . y(t+j −1)   .. ..   h i .  y(t+1)  .   = Γi x(t) . . . x(t+j −1) +   . . ..   ..   y(t+i−1) ... y(t+j +i−2)    u(t) u(t+1) . . . u(t+j −1) e(t)    .. ..    .  u(t+1)   e(t+1) .  + H si  +H i     . .. . .. ..    .    u(t+i−1)

c P. Trnka, 2007 °

...

u(t+j +i−2)

e(t+i−1)

...



..

.

CK

e(t+1) . . . .. . .. . ...

 e(t)

     e(t+1)   . ..     .   I e(t+i−1)

 e(t+j −1)  ..   . ,    e(t+j +i−2)

14

OaF - t´ ema 7


dostáv´ ame tak rovnici rozˇs´ıˇreného stavového modelu, kterou m˚ uˇzeme vyuˇzit´ım Hankelovsk´ ych matic zapsat zkrácenˇe jako Y = Γi X + H i U + H si E

V´ yznam rozˇ s´ıˇ ren´ eho stavov´ eho modelu • spojuje rekurzivn´ı stavov´ y model a experimentáln´ı data do jedn´ e line´ arn´ı rovnice • popisuje vazbu mezi maticemi parametr˚ u (Γi , H i a H si vzniklé z parametr˚ u stavového modelu A, B, C, D, K) a maticemi signál˚ u Y , U , E. • rovnice rozˇs´ıˇreného stavového modelu je v´ ychoz´ım bodem 4SID metod Zaˇcneme-li se na signálové matice d´ıvat jako na ˇrádkové vektory generuj´ıc´ı prostory, potom rozˇs´ıˇren´ y stavov´ y model ukazuje vzájemnou vazbu mezi tˇemito prostory (násoben´ı matic´ı zleva odpov´ıdá lineárn´ı kombinaci ˇrádk˚ u). V´ ystupn´ı prostor row(Y ) tak dostaneme jako lineárn´ı kombinaci prostor˚ u generovan´ ych poˇcáteˇcn´ımi stavy row(X), vstupem row(U ) a ˇsumem row(E). Prostory jsou invariantn´ı ke zmˇenˇe báze stavového modelu. Pro regulárn´ı transformaˇcn´ı matici T ∈ Rn×n a transformaci stav˚ u x0 (t) = T x(t) plat´ı ekvivalence row(X 0 ) = row(T X) = row(X),

col(Γ0i ) = col(Γi T ) = col(Γi ),

row(∆0i ) = row(T ∆i ) = row(∆i ).

c P. Trnka, 2007 °

15

OaF - t´ ema 7

7.5


Deterministick´ a 4SID identifikace

Pro deterministick´ a data staˇc´ı uvaˇzovat rozˇs´ıˇren´ y stavov´ y model bez ˇsum reprezentuj´ıc´ıho ˇclenu H si E Y = Γi X + H i U . Dvˇe hlavn´ı tˇr´ıdy algoritm˚ u 1. projekˇcn´ı 2. pr˚ unikové Liˇs´ı se zp˚ usobem z´ısk´ an´ı prostor˚ u, které jsou dále pouˇzitelné pro odhad parametr˚ u stavového modelu. Ortogon´ aln´ı projekce ve 2D Ortogonáln´ı projekce odpov´ıd´ a rozkladu vektoru do dvou vzájemnˇe kolm´ ych prostor˚ u. Rozklad ˇ r´ adkov´ eho vektoru a ∈ R2 do (prostoru) vektoru b ∈ R2 a ortogonáln´ıho vektoru b⊥ ∈ R2 (tj. pro skal´ arn´ı souˇcin plat´ı b⊥ bT = 0) ⊥

⊥

a = a/b + a/b = l1 b + l2 b

k z´ısk´ an´ı koeficientu l1 vyn´ asob´ıme rovnici zprava bT abT T

ab

c P. Trnka, 2007 °

= l1 bbT + l2 b⊥ bT = l1 bb

b

⊥

a/b ⊥

a

a/b

b

T

16

OaF - t´ ema 7


odtud pro l1 l1 = abT (bbT )−1 a pro ortogonáln´ı projekci a do b a/b = l1 b = abT (bbT )−1 b = aΠb , kde Πb = bT (bbT )−1 b je ortogonáln´ı projekˇcn´ı matice do prostoru generovaného vektorem b. Dosazen´ım l1 do p˚ uvodn´ıho rozkladu dostaneme pro l2 b⊥

l2 b

a =

abT (bbT )−1 b + l2 b⊥

⊥

a[I − bT (bbT )−1 b]

=

a pro ortogonáln´ı projekci a do b⊥ ⊥

a/b = l2 b⊥ = a[I − bT (bbT )−1 b] = aΠb⊥ , kde Πb⊥ = I − bT (bbT )−1 b = I − Πb je ortogonáln´ı projekˇcn´ı matice do prostoru generovaného vektorem bT .

a/b

Shrnut´ı:

⊥

a/b

= aΠb

Πb

=

bT (bbT )−1 b

= aΠb⊥

Πb⊥

=

I − Πb

c P. Trnka, 2007 °

17

OaF - t´ ema 7


Ortogon´ aln´ı projekce ˇ r´ adkov´ ych prostor˚ u matic Uvaˇzujme ˇr´ adkové prostory obecn´ ych matic A ∈ Rm×n , B ∈ Rm×n . Odvozen´ı ortogonáln´ı projekce je zcela shodné s pˇredchoz´ı projekc´ı vektor˚ u ve 2D - ˇrádkové vektory jsou zamˇenˇeny za matice. ˇ adkov´ R´ y prostor matice A lze jednoznaˇcnˇe rozloˇzit do dvou vzájemnˇe kolm´ ych ˇrádkov´ ych prostor˚ u B a B⊥ T ⊥ (plat´ı B B = 0) A = A/B + A/B ⊥ = LB B + LB ⊥ B ⊥ A/B A/B

⊥

=

AΠB

ΠB

=

B T (BB T )−1 B

=

AΠB ⊥

ΠB ⊥

=

I − ΠB

Vˇsimnˇete si, ˇze A/B = AB T (BB T )−1 B = AB + B ⇒ souvislost s u ´lohou LS. Nejprve je ve smyslu LS vyˇreˇsena u ´loha A = θB a následnˇe je projekce dána jako A/B = θB. Pˇ r´ıklad:





A=

1

2

3

3

2

1



Projekˇcn´ı matice:

Projekce matice A:

c P. Trnka, 2007 °



 B=

1 0

0

0 1 0  1  ΠB =   0 0  1 A/B =  3

 0

0



 1 0   0 0  2 0  2 0



ΠB ⊥

A/B

0  =  0 0  ⊥

0

0



 0   0 1  0 0 3  = 0 0 1 0

18

OaF - t´ ema 7

7.5.1


Projekˇ cn´ı deterministick´ y algoritmus

Model Y = Γi X + H i U , ukazuje, ˇze ˇr´ adkov´ y prostor Hankelovy matice v´ ystup˚ u vzniká jako souˇcet prostor˚ u generovan´ ych ˇrádky stavové posloupnosti X pˇres matici pozorovatelnosti Γi a Hankelovy matice vstup˚ u U pˇres matici impulsn´ıch posloupnost´ı H i . Nejjednoduˇsˇs´ı deterministick´ y algoritmus vyuˇz´ıvá skuteˇcnosti, ˇze pro identifikaci parametr˚ u stavového modelu staˇc´ı znalost sloupcového prostoru rozˇs´ıˇrené matice pozorovatelnosti Γi a vstupnˇe/v´ ystupn´ı dat. Matici pozorovatelnosti z´ısk´ ame projekc´ı rovnice rozˇs´ıˇreného stavového modelu do prostoru kolmého na prostor U (pˇren´ asoben´ım projekˇcn´ı matic´ı ΠU ⊥ zprava) Y ΠU ⊥

= (Γi X + H i U ) ΠU ⊥

Y ΠU ⊥

=

Γi XΠU ⊥

t´ım se zbav´ıme ˇclenu H i U , nebot’ U ΠU ⊥ = 0. Z posledn´ı rovnice plyne, ˇze za podm´ınky dostateˇcného vybuzen´ı má matice Y ΠU ⊥ stejn´ y sloupcov´ y prostor jako Γi . Bázi tohoto prostoru z´ıskáme z SVD a ta bude tvoˇrit matici pozorovatelnosti jedné z moˇzn´ ych realizac´ı stavového modelu (nejednoznaˇcnost volby báze). Singulárn´ım rozkladem   h i S1 0 h iT ³ ´³ ´ 1/2 1/2 T  V1 V2 Y ΠU ⊥ = U SV = U 1 U 2  = U 1S1 S 1 V T1 0 0

c P. Trnka, 2007 °

19

OaF - t´ ema 7


kde matice U , S a V T jsou rozdˇeleny podle poˇctu nenulov´ ych singulárn´ıch ˇc´ısel. Dostaneme 1/2

Γi = U 1 S 1

Poˇcet nenulov´ ych singulárn´ıch ˇc´ısel odpov´ıdá ˇrádu systému n. Urˇ cen´ı parametr˚ u A, B, C, D ze znalosti Γi a vstupnˇ e/v´ ystupn´ıch dat Stejnˇe jako v realizaˇcn´ım algoritmu lze A a C urˇcit z invariantnosti Γi v˚ uˇci posunu +

A

=

Γi Γi ,

C

=

Γi (1 : n, :).

+ Urˇcen´ı C a D nen´ı tak pˇr´ımoˇcaré. Pˇren´ asoben´ım rozˇs´ıˇreného stavového modelu Γ⊥ zprava dostaneme i zleva a U + Γ⊥ = Γ⊥ i YU i H i.

Oznaˇc´ıme k-t´ y sloupec levé strany rovnice jako Mk a k-t´ y sloupec Γi jako Lk . Rovnici pak m˚ uˇzeme pˇrepsat na   D    h i h i D  CB   M1 . . . Mi = L1 . . . Li  ..   ..   . . CB   i−2 CA B . . . CB D

c P. Trnka, 2007 °

20

OaF - t´ ema 7


Nasklád´ an´ım sloupc˚ u levé strany rovnice nad sebe dostaneme     L1 L2 . . . Li     M1 .. .    L  . ..  I 0 D 2  ..     ,  . =     .. B  ...  0 Γi .   Mi Li 0 coˇz je lineárn´ı soustava rovnic v promˇenn´ ych D a B   D , M = L B jej´ımˇz ˇreˇsen´ım dostaneme jediné neznámé B a D, nebot’ M a L jsou funkcemi I/O dat a Γi .

c P. Trnka, 2007 °

OaF - t´ ema 7

7.5.2

21


Pr˚ unikov´ y deterministick´ y algoritmus

Zaloˇzen´ y na moˇznosti z´ıskat posloupnost stav˚ u z pr˚ uniku dvou prostor˚ u, které jsou generovány I/O daty. Rozdˇelen´ı signálov´ ych matic U , Y a E na ˇcást historickou (past) a budouc´ı (future)   u(0) u(1) ... u(j − 1)   i . . . délka okna do historie   u(1) u(2) ... u(j)     h . . . délka okna do budoucnosti . . . .   .. .. .. ..   j = N − i + 1 kde N je poˇcet vzork˚ u       Up u(i − 1) u(i) . . . u(i + j − 2)   =    i ¿ j, h À j ⇒ ˇsiroké obdéln´ıkové matice Uf u(i) u(i + 1) . . . u(i + j − 1)       u(i + 1)  u(i + 2) . . . u(i + j)   i a h se vol´ı o nˇeco vˇetˇs´ı neˇz oˇcekávan´ y   . . . .   . . .. .. . .   ˇrád systému u(i + h − 1) u(i + h) . . . u(i + h + j − 2) podobnˇe dostaneme Y p , Y f , E p a E f . Rozdˇelen´ı stavové posloupnosti h i h i X p = x(0) . . . x(j − 1) , X f = x(i) . . . x(j + i − 1) Stavovou posloupnost X f pak dostaneme jako pr˚ unik historick´ ych a budouc´ıch dat     Uf Up T        row row (X f ) = row Yp Yf To ukazuje, ˇze stav pˇredstavuje propojen´ı mezi historick´ ym a budouc´ım v´ yvojem systému. c P. Trnka, 2007 °

22

OaF - t´ ema 7


Nalezen´ı pr˚ uniku vych´ az´ı z tzv. principi´ aln´ıch u ´ hl˚ u a principi´ aln´ıch smˇ er˚ u mezi lineárn´ımi prostory, coˇz je zobecnˇen´ı koncepce u ´hlu mezi dvˇema vektory na u ´hly mezi dvˇema prostory. Pˇ r´ıklad: dvˇe roviny A a B ve 3D prostoru (obrázek) Prvn´ı principiáln´ı smˇery hledáme jako jednotkové vektory a1 ∈ A a b1 ∈ B, které sv´ıraj´ı nejmenˇs´ı u ´hel. Jsou to vektory leˇz´ıc´ı v prostoru pr˚ uniku a sv´ıraj´ı prvn´ı principiáln´ı u ´hel θ1 = 0. Druhé principiáln´ı smˇery hledáme uˇz v menˇs´ım prostoru a2 ∈ A − ha1 i a b2 ∈ B − hb1 i a opˇet pod nejmenˇs´ım u ´hlem θ2 = θ. Prostor pr˚ uniku je potom urˇcen principi´ aln´ımi smˇ ery odpov´ıdaj´ıc´ı nulov´ ym principi´ aln´ım u ´ hl˚ um.

A a2 a1=b1 θ1=0

θ2

B

b2

Praktické nalezen´ı principiáln´ıch u ´hl˚ u a smˇer˚ u a pr˚ uniku umoˇzn ˇuje SVD. Oznaˇc´ıme W p jako historická data a W f jako budouc´ı data     Up Uf   Wp =  Wf =  Yp Yf Singulárn´ım rozkladem rozloˇz´ıme souˇcin ortogonáln´ıch projekˇcn´ıch matic   h i S1 0 h iT T  V1 V2 ΠWp ΠWf = U SV = U 1 U 2  0 0 c P. Trnka, 2007 °

23

OaF - t´ ema 7


Bázi prostoru pr˚ uniku pak dostaneme z lev´ ych (stejnˇe tak i z prav´ ych) singulárn´ıch vektor˚ u pˇr´ısluˇsn´ ych jedniˇckov´ ym singulárn´ım ˇc´ısl˚ um (singulárn´ı ˇc´ısla odpov´ıdaj´ı kos´ın˚ um principiáln´ıch u ´hl˚ u) X f = U T1 . Urˇ cen´ı parametr˚ u A, B, C, D ze znalosti posloupnosti stav˚ u a I/O dat N −1

Známe-li posloupnosti {u(t), x(t), y(t)}0 potom stavov´ y model dává pˇr´ımo vazbu mezi tˇemito daty      x(0) . . . x(N − 2) x(1) . . . x(N − 1) A B .   = u(0) . . . u(N − 2) C D y(0) . . . y(N − 2) Dostáv´ ame tak pˇreurˇcenou soustavu line´ arn´ıch rovnic, ze které m˚ uˇzeme urˇcit A, B, C a D. Dostateˇ cn´ e vybuzen´ı (Persistency of excitation) N −1 Posloupnost {u(t)}0 je dostateˇcnˇe vybuzuj´ıc´ı ˇrádu ` jestliˇze jej´ı Hankelova matice   u(0) u(1) . . . u(N − ` − 2)   ..    u(1)  .     .. ..   . .   u(` − 1) . . . . . . u(N − 1) má hodnost rovnu `. (Neplat´ı napˇr´ıklad pro periodické signály)

c P. Trnka, 2007 °

24

OaF - t´ ema 7

7.6


Stochastick´ a 4SID identifikace

V pˇr´ıtomnosti ˇsumu bude stavov´ a rovnice rozˇs´ıˇreného stavového modelu Y = Γi X + H i U + H si E. Jak se zachovaj´ı pˇredchoz´ı algoritmy: • U projekˇcn´ıho algoritmu bude s nar˚ ustaj´ıc´ım ˇsumem r˚ ust poˇcet nenulov´ ych singulárn´ıch ˇc´ısel Y ΠU ⊥ (to zt´ıˇz´ı aˇz znemoˇzn´ı urˇcen´ı ˇr´ adu) a odhad bude vych´ ylen´ y i pro nekoneˇcnˇe velkou mnoˇzinu I/O dat, nebot’ projekc´ı rozˇs´ıˇreného stavového modelu do U ⊥ Y ΠU ⊥

=

Γi XΠU ⊥ + H i U ΠU ⊥ + H si EΠU ⊥

Y ΠU ⊥

=

Γi XΠU ⊥ +

+H si EΠU ⊥

neodstran´ıme ˇsumem vznikl´ y ˇclen H si E, kter´ y je ortogonáln´ı k U . • U pr˚ unikového algoritmu se bude ztrácet pr˚ unik mezi ˇrádkov´ ymi prostory W p a W f - pˇresto principiáln´ı smˇery pˇr´ısluˇsné nejmenˇs´ım u ´hl˚ um dávaj´ı dobrou aproximaci stav˚ u pro rostouc´ı mnoˇzinu I/O dat.

c P. Trnka, 2007 °

25

OaF - t´ ema 7

7.6.1


N4SID (Numerical 4SID)

Algoritmus zaloˇzen´ y na ˇ sikm´ e projekci umoˇzn´ı z rovnice rozˇs´ıˇreného stavového modelu asymptoticky odstranit jak vliv deterministického vstupu, tak vliv ˇsumu. ˇ Sikm´ a projekce ˇ Obecné matice A ∈ Rp×j , B ∈ Rq×j a C ∈ Rr×j . Sikm´ a projekce rozklád´ a matici A do tˇr´ı ˇc´ ast´ı ⊥  B  A = A / C + A / B + A/  B C C a je definov´ ana jako h A/C =A C

T

B

T

i

 

B

+    h i I  C T B T   r×r  C 0 B C

A A/ B

B

A/

C

A/C B

B C

C

nebo h ih i+ A / C = A/B ⊥ C/B ⊥ C. B

ˇ V´ yznam: Sikm´ a projekce pˇredpoklád´ a, ˇze A vzniklo jako lineárn´ı kombinace C, B a jeˇstˇe neznámé sloˇzky (ˇsum). V´ ysledkem projekce A / C je jak´ y byl na této lineárn´ı kombinaci pod´ıl sloˇzky C. B c P. Trnka, 2007 °

26

OaF - t´ ema 7


ˇ Sikmou projekc´ı rozˇs´ıˇreného stavového modelu do W p podél U f dostaneme Y f / Wp Uf

= (Γh X f + H h U f + H sh E f ) / W p Uf

= Γh X f / W p + H h U f / W p + H sh E f / W p . Uf

Uf

Uf

pro jednotlivé ˇcleny: • Pro poˇcet dat a délku okna do historie jdouc´ı k nekoneˇcnu (i, j → ∞) dostáváme projekc´ı X f do W p pˇres ˆ f z historick´ U f nejlepˇs´ı odhad stav˚ uX ych dat ˆf. lim X f / W p = lim X f /W p = X

i,j→∞

Uf

i,j→∞

• Z definice ˇsikmé projekce U f / W p = 0. Uf

• Pro nezávislé ˇsumy a pro poˇcet dat jdouc´ı k nekoneˇcnu lim E f / W p = 0.

j→∞

Uf

c P. Trnka, 2007 °

27

OaF - t´ ema 7


Shrnut´ı: ˆf lim Y f / W p = Γh X

i,j→∞

Uf

Ze ˇsikmé projekce tak m˚ uˇzeme singulárn´ım rozkladem odhadnout ˇrád systému n, rozˇs´ıˇrenou matici pozorovatelnosti Γh a stavovou posloupnost X f . Algoritmus: 1. SVD rozklad ˇsikmé projekce Y f / W p = U SV T = Uf

h

i U1

U2

 

 S1

0

0

S2



h

iT V1

V2

2. Inspekc´ı singulárn´ıch ˇc´ısel v S odhadneme ˇrád systému n a dle nˇej rozdˇel´ıme U a V T . Dostaneme odhady 1/2

Γh = U 1 S 1 ,

ˆ f = S 1/2 V T1 X 1

3. Z odhadnuté stavové posloupnosti a I/O dat odhadneme parametry stavové modelu A, B, C, D stejnˇe jako v projekˇcn´ım algoritmu      A B x ˆ(0) . . . x ˆ(N − 2) x ˆ(1) . . . x ˆ(N − 1)  + ε.   = u(0) . . . u(N − 2) C D y(0) . . . y(N − 2) 4. Kalmanovo zes´ılen´ı K a kovarianci ˇsumu Re odhadneme z rezidu´ı ε   Σ11 Σ12 K = Σ12 Σ−1 1 22 , = εεT kde Σ =  j − (n + m)(n + l) Σ21 Σ22 Re = Σ22 ,

c P. Trnka, 2007 °

28

Subspace Identification Methods

Recommend Documents