Akademia Górniczo-Hutnicza
Katedra Robotyki i Mechatroniki
Identyfikacja i analiza sygnałów
Laboratorium 9
Metody błędów predykcji
Modelowanie z zastosowaniem metody błędów predykcji
Metoda najmniejszych kwadratów, której była mowa na poprzednim laboratorium, jest z pew-
nością prosta w użyciu, chociaż uzyskanie estymatorów zgodnych wymaga przyjęcia stosun-
kowo silnych założeń. W pewnych sytuacjach brak zgodności może być jednak tolerowany. Tak
na przykład jeśli stosunek sygnału do szumu jest wysoki, to obciążenie będzie małe. Przy pro-
jektowaniu regulatorów dopuszczalne jest pewne obciążenie estymatora modelu, gdyż popraw-
nie zaprojektowany układ regulacji posiada pewną odporność na zmiany parametrów modelu.
Są jednak sytuacje, w których estymatory zgodne mają wielkie znaczenie. Metodę najmniej-
szych kwadratów można zmodyfikować tak aby dawała estymatory zgodne. Są to następujące
modyfikacje:
1. Minimalizacja błędów predykcji dla bardziej złożonych struktur modeli parametrycz-
nych. Prowadzi to do klasy metod błędu predykcji, którą będziemy się zajmować w tym
rozdziale.
2. Modyfikacja równań normalnych, odpowiadających estymatorowi metody najmniej-
szych kwadratów. To z kolei prowadzi do metod zmiennej instrumentalnej.
W tym laboratorium zajmiemy się metodą błędów predykcji. Jej schemat można opisać w na-
stępujący sposób.
1. Wybierz strukturÄ™ modelu w postaci:
y(t) =ð G(q-ð1;qð)u(t) +ð (q-ð1;qð)e(t)
i predykator w postaci
w
(t | t -ð1;qð) =ð L1(q-ð1;qð)y(t) +ð L2(q-ð1;qð)u(t)
który zależy wyÅ‚Ä…cznie od przeszÅ‚ych danych tylko wtedy, filtry L1(0;qð)=0 i L2(0;qð)=0
2. Wybierz funkcjÄ™ celu h(Q)
3. Wyznacz ocenÄ™ parametrów qð ðjako globalny punkt minimalizujÄ…cy funkcjÄ™ strat h(RN(qð))
Schemat metody błędów predykcji pokazano na rysunku poniżej:
Modele regresyjne wykorzystują najczęściej operatorowy zapis dynamiki systemu w postaci
transmitancji (funkcji przejścia) lub równań stanu. Przyjęto zapis obiektu w postaci torów
wejść-wyjść oraz zakłóceń-wyjść. Oznacza to, że dla modeli wielomianowych zakłócenia opi-
suje się jako sprowadzone na wyjście obiektu (Rys.1).
y(t) =ð G(s)u(t) +ð (s)e(t) (1)
lub
y(i) =ð G(z-ð1)u(i) +ð (z-ð1)e(i) (2)
Jeżeli zakłócenia oddziałujące na obiekt nie są modelem białego szumu, przyjmuje się opis zło-
żony w postaci filtru H(z-1) na którego wejściu jest biały szum. Dyskretne transmitancje odpo-
wiadające torom związanym z modelowaniem zależności pomiędzy wejściem i wyjściem G(z-1)
oraz zależności pomiędzy zakłóceniami i wyjściem H(z-1) można wyrazić jako funkcje wymier-
ne operatora z-1 (dla modeli z czasem ciągłym operatora s).
Rys. 1. Tory wejścia-wyjścia (tor sterowania) oraz zakłócenia-wyjścia (tor
zakłóceń) (modeli obiektów w ciągłej i dyskretnej dziedzinie czasu
Modele obiektów mogą być definiowane za pomocą następujących wielomianów operatoro-
wych:
A(z-ð1) =ð1+ð a1z-ð1 +ð a2z-ð2 +ð, ...,+ðanAz-ðnA, (3)
B(z-ð1) =ð b0 +ð b1z-ð1 +ð b2z-ð2 +ð, ...,+ðbnBz-ðnB , (4)
C(z-ð1) =ð1+ð c1z-ð1 +ð c2z-ð2 +ð, ...,+ðcnCz-ðnC , (5)
D(z-ð1) =ð1+ð d1z-ð1 +ð d2z-ð2 +ð, ...,+ðdnDz-ðnD, (6)
F(z-ð1) =ð1+ð f1z-ð1 +ð f2z-ð2 +ð, ...,+ð fnF z-ðnF . (7)
Dodatkowo wprowadza się modele zakłóceń deterministycznych (trendy wielomianów
dowolnego stopnia, np. liniowy, kwadratowy, sześcienny, trendy okresowe lub stała ) oraz
zakłócenia stochastyczne (dla C(z-1)=1 i m=1, szczególny przypadek niestacjonarności tzw.
błądzenia przypadkowego) w postaci
1 dla i =ð 0
d ìð
d(i) =ð dð (i), dð (i) =ð (8)
íð0 dla i Ä…ð 0.
Ag (z-ð1)
îð
m
Ńðm =ð (ð1-ð z-ð1)ð . (9)
gdzie wielomian Ag(z-1) generujący zakłócenia może być przyjęty przykładowo w postaci
Ag (z-ð1) =ð1-ð 2z-ð1 cos(wðs ) +ð z-ð2 . (10)
co odpowiada zależnoÅ›ci czasowej d(i), gdzie wðs jest czÄ™stoÅ›ciÄ… sinusoidy zakłócajÄ…cej
d(i) =ð d sin(wðsi) . (11)
Poprzez niestacjonarność rozumie się dryft wybranych wartości zmiennych w czasie (zmiana
średniej, wariancji, momentów wyższych rzędów). Można wyróżnić następujące jej rodzaje:
a) niestacjonarność deterministyczną w postaci nałożonego trendu wielomianowego (8),
b) niestacjonarność stochastyczną związaną z pierwiastkami jednostkowymi (9)
Uwzględniając wpływ wszystkich możliwych zakłóceń, uogólniony model przedstawia
schemat znajdujÄ…cy siÄ™ na Rys. 2.
Rys. 2. Uogólniony model dynamiczny obiektu
Wielomianowy model ogólny SISO zapisujemy następująco po podstawieniu do (2)
poszczególnych wielomianów
B(z-ð1) C(z-ð1)
A(z-ð1)[ðŃðm y(z-ð1)]ð=ð [ðŃðmu(z-ð1)]ð+ð e(z-ð1) (12)
F(z-ð1) D(z-ð1)
W celu budowy modelu systemu dynamicznego należy określić przyszłą strukturę modelu
(Tabela 1). Dokonujemy tego poprzez włączenie wybranych kanałów (transmitancji operatora
z-1) modyfikujących własności ciągów czasowych na wejściu modelu.
Przyjmuje się, że skalarny dyskretny identyfikowalny model obiektu dynamicznego opisuje
następująca struktura
Tabela 1. Zestawienie wybranych liniowych modeli regresyjnych
ARMAX(nA, nB, nC, d, m) (13)
Istotnym elementem modelu obiektu fizycznego jest dobranie odpowiedniego opóz
nienia
pomiędzy wejściem i wyjściem modelu. Przykładowo dla obiektów mechanicznych należy
uwzględnić prędkość rozchodzenia się odkształceń ośrodka (m. in. wymuszanie drgań falą
akustyczną). Orientacyjnie można wyznaczyć krotność opóz
nienia wynikającego z przyjętej
częstotliwości próbkowania zaokrąglając wartość następującego wyrażenia
l
d =ð (14)
v ×ðTs
gdzie: d jest szacowanym opóz
nieniem, l odległość w ośrodku fizycznym pomiędzy punktem
przyłożenia wymuszenia oraz punktem pomiaru odpowiedzi obiektu, v jest prędkością
rozchodzenia się odkształceń ośrodka, Ts jest okresem próbkowania.
Metody identyfikacji
Stosowane są dwie podstawowe klasy metod identyfikacji parametrów modelu. Metoda
błędów równań wielomianów (omawiana na poprzednim laboratorium), w której
identyfikowane parametry są zależne liniowo od błędu estymacji oraz metoda błędu wyjścia z
nieliniowymi zależnościami parametrów od błędów predykcji. Metoda pierwsza stosowana
jest do identyfikacji wielomianowych struktur ARX natomiast metoda druga do identyfikacji
modelu ogólnego PEM (ang. Predict Error Method) oraz modeli pochodnych wynikających z
uproszczenia modelu PEM, są to m. in. modele BJ, OE, ARMAX, jak również modele w
innowacyjnej formie równań stanu.
Rys. 3. Estymacja parametrów modelu w przypadku metody błędu równań wielomianów
Rys. 4. Estymacja parametrów modelu w przypadku metody błędu wyjściowego
Rys. 3 i Rys. 4 przedstawiają schematycznie obydwie wymienione metody. W metodzie błędu
równań wielomianów w celu uzyskania estymatora o pożądanych właściwościach (zgodny,
zbieżny i efektywny) wymaga się niewielkich zakłóceń oddziałujących na model.
Przedstawiono zakłócenia oddziałujące na wejście ein oraz na wyjście eout do których można
zaliczyć modele szumów pomiarowych (miedzy innymi błędy grube itp.). Innym typem
zakłóceń są oznaczone jako e błędy związane z oddziaływaniem szumów na zmienne stanu, o
ile jest przyjęty taki model właściwy obiektu (ang. true system).
Modele torów zakłóceń dla obiektów mechanicznych
Można rozważać obiekt mechaniczny o parametrach fizycznych rozłożonych przestrzenie.
Wiadomo, że pobudzenie takiego obiektu do drgań w punkcie p1 oraz p2 powoduje
odpowiedz
, którą można obserwować w wielu punktach obiektu. Odpowiedz
w punkcie p3
stanowi wypadkową zależną od wymuszeń u1 i u2.
Rys. 5. Schemat systemu rzeczywistego z naniesionymi punktami, w któ-
rych przyłożono wymuszenie oraz mierzono odpowiedz

W wielu praktycznych przypadkach identyfikuje siÄ™ uproszone liniowe modele SISO
rzeczywistych obiektów mechanicznych, przyjmując jako sygnał wejściowy i wyjściowy
wybrane wielkości fizyczne dostępne pomiarowo.
Rys. 6. Struktura modelu regresyjnego SISO, wynikajÄ…ca z uproszczenia
modelu MIMO
Zastosowanie modelu ARX oraz OE wymaga przyjęcia odpowiednich stopni wielomianów
A(z-1) dla modelu ARX oraz F(z-1) dla modelu OE pozwalających opisać składowe
harmoniczne widma. Z punktu widzenia interpretacji fizycznej jest to niepoprawne ponieważ
zakłada się, że drgania w punkcie y1 powodowane są wyłącznie wymuszeniem u1. W takich
przypadkach modele BJ oraz PEM o bardziej złożonym sposobie modelowania zakłóceń dają
większe możliwości parametryzacji umożliwiające uzyskanie estymatora zgodnego, poprzez
separację zakłóceń od toru wejście-wyjście (toru sterowania). Umożliwiają również
zachowanie interpretacji fizycznej właściwości obiektu. W przypadku obiektu z Rys.
uwzględniając wiedzę a priori wiadomym jest, że na odpowiedz
w punkcie y1 wpływa
niemierzalne wymuszenie u2 można ten fakt uwzględnić wprowadzając wspólny wielomian
A(z-1) dla modelu PEM interpretowany jako część wspólna mianownika toru wejście-wyjście
oraz toru zakłócenia-wyjście co prowadzi do modelu o strukturze PEM(nA, nB, nC=1, nD,
nF=1,k).
Rys. 7. Szczególny przypadek struktury PEM
Wykazano, że każdy obiekt można zidentyfikować przy użyciu najprostszej struktury ARX,
jednak o większej liczbie parametrów niż struktury złożonych modeli.
Szczególne zalety modelu BJ w identyfikacji obiektów mechanicznych
Model o strukturze BJ umożliwia rozprzęgniecie mianowników funkcji przejścia dla sygnału
użytecznego G(z-1) oraz dla szumu H(z-1). Przykład odmiennych właściwości funkcji
przejścia kształtujących sygnał użyteczny oraz szum przedstawia Rys. 8 zawierających
charakterystyki widmowe wykonane na podstawie modelu BJ(12,8,8,12,1) dla którego
stopnie poszczególnych wielomianów wynoszą: nB=12, nC=8, nD=8, nF=12 oraz opóz
nienie
k=1. Rozprzęgnięcie mianowników ma duże znaczenie podczas identyfikacji obiektów
charakteryzujących się nieliniowościami, wtedy ustalając nF << nD można opisać zakłócenia
charakteryzujące się dużą liczbą harmonicznych częstotliwości własnych.
Rys. 8. Funkcje przejścia dla poszczególnych wielomianów modelu związanych z: a)
sygnałem użytecznym G(z-1) oraz b) z zakłóceniami H(z-1).
PRZYKA
AD
Identyfikowano obiekt o transmitancji
1 1
H (s) =ð +ð
s2 +ð10s +ð (2pð 30)2 s2 +ð 2s +ð (2pð 80)2
którego wyjście przekształca nieliniowa funkcja statyczna fN
fN (ð.)ð=ð fN[ð109 ×ð y(t)3]ð
éð Å‚ð
2s2 +ð12s +ð 288200
y(s) =ð fN Ä™ðs4 u(s)Å›ð
+ð12s3 +ð 288200 s2 +ð 2598000 s +ð 8.977 ×ð108 ûð
ëð
Obiekt pobudzany jest sygnałem harmonicznym o modelu
u(t) =ð sin(2pð 60t)
gdzie tÎð(0,0.8) [s], czÄ™stotliwość próbkowania fp=500 [Hz]
éð Å‚ð
3.798 ×ð10-ð6 z3 -ð1.695 ×ð10-ð6 z2 -ð1.705 ×ð10-ð6 z +ð 3.738 ×ð10-ð6
y(z) =ð fN Ä™ð u(z)Å›ð
z4 -ð 2.911z3 +ð 3.945z2 -ð 2.882z +ð 0.9763
ëð ûð
Rys. 3. Przebiegi czasowe wejścia-wyjścia symulowanego obiektu dla przypadku
liniowego (linia ciągła) oraz nieliniowego (linia kreskowa)
Rys. 4. Widma mocy w skali logarytmicznej przebiegów czasowych wejścia-wyjścia
symulowanego obiektu dla przypadku liniowego (linia ciągła) oraz nieliniowego (linia
kreskowa)
Rys. 5. Przykłady identyfikacji symulowanego obiektu nieliniowego (23) modelem
ARX(4,3,1) (lewa strona) oraz modelem ARX(6,3,1) (prawa strona). Linia ciągła oznacza
model dyskretny (23), linia kreskowa model ARX, linia kreska-kopka przedział ufności na
poziomie 1-að=0.95, amplituda w skali logarytmicznej
Rys. 6. Przykłady identyfikacji symulowanego obiektu nieliniowego (23) modelem
ARX(4,3,1) (lewa strona) oraz modelem ARX(6,3,1) (prawa strona) dla toru zakłóceń H(z-1).
Linia ciągła oznacza model dyskretny (23), Linia kreska-kopka przedział ufności na poziomie
1-að=0.95, amplituda w skali logarytmicznej.
Rys. 7. Przykłady identyfikacji symulowanego obiektu nieliniowego (23) modelem
ARMAX(4,3,3,1) (lewa strona) oraz modelem OE(3,4,1) (prawa strona). Linia ciągła oznacza
model dyskretny (23), linia kreskowa model ARMAX/OE, linia kreska-kopka przedział
ufnoÅ›ci na poziomie 1-að=0.95, amplituda w skali liniowej
Rys. 8. Przykłady identyfikacji symulowanego obiektu nieliniowego (23) modelem
ARMAX(4,3,3,1) (lewa strona) oraz modelem OE(3,4,1) (prawa strona) dla toru zakłóceń
H(z-1). Linia kreska-kopka przedziaÅ‚ ufnoÅ›ci na poziomie 1-að=0.95, amplituda w skali
logarytmicznej.
Rys. 9. Przykłady identyfikacji symulowanego obiektu nieliniowego (23) modelem
BJ(3,2,8,4,1) (lewa strona) oraz modelem BJ(4,2,8,4,1) (prawa strona). Linia ciągła oznacza
model dyskretny (23), linia kreskowa model BJ, linia kreska-kopka przedział ufności na
poziomie 1-að=0.95, amplituda w skali liniowej
Rys. 10. Przykłady identyfikacji symulowanego obiektu nieliniowego (23) modelem
BJ(3,2,8,4,1) oraz modelem BJ(4,2,8,4,1) dla toru zakłóceń H(z-1) (identyczne stopnie
wielomianów toru zakłócenia). Linia kreska-kopka przedziaÅ‚ ufnoÅ›ci na poziomie 1-að=0.95,
amplituda w skali logarytmicznej.
Rys. 11. Przykłady symulacji wyjścia na podstawie zidentyfikowanego modelu obiektu
nieliniowego (23) o strukturze BJ(3,2,8,4,1). Linia kreska-kopka przedział ufności na
poziomie 1-að=0.95, amplituda w skali logarytmicznej.
Rys. 12. Przykłady błędu symulacji wyjścia na podstawie zidentyfikowanego modelu
obiektu nieliniowego (23) o strukturze BJ(3,2,8,4,1). Linia kreska-kopka przedział ufności na
poziomie 1-að=0.95, amplituda w skali logarytmicznej.
Wnioski i podsumowanie
Zakładając możliwość identyfikacji modelu SISO uproszczonego w stosunku do rozważanego
modelu MIMO istnieją następujące możliwości parametryzacji torów zakłóceń:
a) model ARX opisujący zakłócenia jako proporcjonalne do wyjścia sygnału użyteczne-
go,
b) model ARMAX wprowadzający częściowo niezależny opis zakłóceń w postaci zer
licznika C(z-1) w torach zakłócenia-wyjścia, mianownik jest taki sam jak w torach
wejścia-wyjścia,
c) model OE pomijający opis zakłóceń przy pomocy funkcji przejścia (na wyjściu tylko
biaÅ‚y szum o wariancji lð),
d) model BJ wprowadzający całkowicie niezależny opis zakłóceń od torów wejścia-
wyjścia w postaci zer licznika C(z-1) oraz biegunów mianownika D(z-1) w torach za-
kłócenia-wyjścia,
e) model PEM wprowadzający całkowicie niezależny opis zakłóceń od torów wejścia-
wyjścia jak w modelu BJ oraz dodatkowo wspólną część w postaci biegunów mia-
nownika A(z-1) dla torów wejścia-wyjścia oraz zakłócenia-wyjścia.
W czasie działania obiektu przemysłowego można wyróżnić dwie klasy modelowanych
zmian:
a) zmiany otoczenia w którym działa obiekt,
b) zmiany właściwości obiektu
Z praktycznego punktu widzenia prowadzi to do błędów stosowania modelu obiektu, jeżeli w
otoczeniu pojawi siÄ™ z
ródło dodatkowych zakłóceń lub w obiekcie wystąpią zmiany
właściwości. W celu zapewnienia zgodności estymatora należy poprawnie dobrać strukturę
modelu wykorzystując przede wszystkim wiedzę a priori dotyczącą zakłóceń otoczenia oraz
właściwości obiektu. W wielu przypadkach najprostsza struktura ARX nie zapewnia nawet
dostatecznej dokładności w odtwarzaniu właściwości obiektu oraz zakłóceń otoczenia.
Dobór struktury w przypadku modeli bardziej złożonych: ARMAX, BJ, PEM jest bardzo
trudnym zadaniem. Niewielkie zmiany w stopniach wielomianów A,B,C,D,F powodują duże
zmiany jakościowe w wynikach stosowania modelu.
W przypadku diagnostycznych zastosowań modelu można wykorzystać zmiany
strukturalne do diagnostyki, np. separacji zmian w nieliniowości obiektu, pojawiania się
ubocznych skutków działania obiektu itd. Wprowadzenie modeli BJ (również PEM) z
oddzielnymi biegunami w torach wejścia-wyjścia oraz zakłócenia-wyjście predysponuje te
struktury do identyfikacji obiektów analizowanych pod kątem drgań.
Zadania do samodzielnego wykonania:
1. Zamodeluj w Simulinku układ z rysunku poniżej i przeprowadz
jego symulacjÄ™ wymusza-
jąc go szumem losowym. Przebiegi czasowe zarówno wejścia (siła f1) jak i wyjścia
(przemieszczenia x2) zapisz do przestrzeni roboczej Matlaba.
Dane: M1 = 10, M2 = 7, C1 = 0.7, C2 = 0.9, K1 = 120, K2 = 200
2. Dla zapisanych danych dokonaj identyfikacji parametrów modeli klasy ARMAX, BJ, OE
i PEM przyjmując porównywalne rzędy odpowiednich wielomianów i zerowe opóz
nienie.
3. Dokonaj symulacji otrzymanych modeli za każdym razem stosując ten sam sygnał wymu-
szający co w zadaniu 1 (siła f1). Porównaj otrzymane wyniki wyliczając współczynnik ko-
relacji pomiędzy przebiegami z kolejnych symulacji modeli parametrycznych, a modelem
referencyjnym z Simulinka.
4. Powtórz zadania 2 i 3 dodając do danych symulacyjnych szum o amplitudzie 0.1 sygnału.