Akademia Górniczo-Hutnicza
Katedra Robotyki i Mechatroniki
Identyfikacja i analiza sygnałów
Laboratorium 11
Ocena jakości estymatorów
Ocena jakości estymatorów w identyfikacji obiektów dynamicznych
Przy identyfikacji modeli układów dynamicznych, oprócz znalezienia (wyestymowania) para-
metrów modelu, musimy również sprawdzić na ile model przez nas uzyskany jest poprawny.
Znaczy to jak dobrze naśladuje rzeczywisty układ, który modelujemy. Służy do tego szereg na-
rzędzi, które są stosowane w zależności od typu modelu. W czasie niniejszego laboratorium stu-
denci zapoznają się z następującymi sposobami oceny jakości estymatorów:
- dla modeli parametrycznych: kryterium Akaike, funkcja kowariancji i wyznaczanie przedzia-
łów ufności,
- dla modeli nieparametrycznych: funkcja koherencji i wskaz
niki FRAC i RVAC,
- dla modeli modalnych: kryterium MAC i pochodne CoMAC i MACCo.
Modele parametryczne
Kryterium informacyjne Akaike  zaproponowane przez Hirotsugu Akaike pod nazwÄ…
kryterium informacyjne w 1971 roku. Jest miarÄ… dobroci dopasowania modeli statystycznych.
Bazuje na pojęciu entropii w efekcie pokazujące względną miarę utraty informacji, gdy dany
model jest wykorzystywany do opisu rzeczywistości obrazuje związek pomiędzy stopniem
złożoności modelu a jego dokładnością. AIC kryterium jest testem pozwalającym na wybór
najlepszego modelu. Dla danego zestawu danych tworzymy grupę modeli, które oceniamy
przy pomocy AIC i wybieramy najlepszy czyli ten, który ma najniższą wartość AIC. Akaike
zaproponował, aby wybierać ten model dla którego najmniejsza jest wartość:
AIC = - 2k - 2ln(L)
gdzie:
" k - liczba parametrów modelu.
" L  estymowane prawdopodobieństwo, przy założeniach danego modelu, uzyskania
takiej właśnie wartości obserwacji j jaka była naprawdę uzyskana
PRZYKA
AD
Dla układu o jednym stopniu swobody i parametrach m,c,k równych odpowiednio 7, 0.7, 700
m=7, c=0.7, k=700
zbuduj model w przestrzeni stanów,
A=[0 1;-k/m -c/m];
B=[0; 1/m];
C=[1 0];
D=0;
SSmodel=ss(A,B,C,D);
dokonaj jego symulacji przy wymuszeniu losowym,
t=[0:0.01:5];
u=3*rand(size(t));
y=lsim(SSmodel, u, t);
otrzymane wyniki zaszumiamy szumem losowym o rozkładzie normalnym, średniej zero i
amplitudzie równej 0.1 amplitudy zakłócanego sygnału,
for c=1:length(u)
u(c)=u(c)*(0.9+0.1*randn(1));
end;
for c=1:length(y)
y(c)=y(c)*(0.9+0.1*randn(1));
end;
dla tak przygotowanych sygnałów dokonujemy identyfikacji modeli typu ARX, ARMAX, BJ
i PEM,
Dane=iddata(y,u',0.01);
ARXmodel=arx(Dane,[2,2,0]);
ARMAXmodel=armax(Dane,[2,2,2,0]);
BJmodel=bj(Dane,[2,2,2,2,0]);
PEMmodel=pem(Dane,[2,2,2,2,2,0]);
następnie przy pomocy kryterium AIC sprawdzamy, który model jest najlepszy:
KryteriumAIC=aic(ARXmodel,ARMAXmodel,BJmodel,PEMmodel)
Powyższy kod pozwala na uzyskanie następującego wyniku:
KryteriumAIC =
-15.4129 -16.2133 -15.8715 -15.8103
Wynika z tego, iż dla tak prostego modelu i przyjętych rzędów modeli najlepszym modelem
jest ARMAX.
Innym sposobem oceny modeli parametrycznych jest wyznaczenie dla nich przedziałów ufno-
Å›ci. Niech zmienna X ma rozkÅ‚ad w populacji z nieznanym parametrem ¸. Z populacji wybie-
ramy próbÄ™ losowÄ… (X , X , ..., X ). PrzedziaÅ‚em ufnoÅ›ci (¸ - ¸ , ¸ + ¸ ) o współczynniku uf-
1 2 n 1 2
noÅ›ci 1 - Ä… nazywamy taki przedziaÅ‚ (¸ - ¸ , ¸ + ¸ ), który speÅ‚nia warunek:
1 2
P(¸ < ¸ < ¸ ) = 1 - Ä…
1 2
gdzie: ¸ i ¸ sÄ… funkcjami wyznaczonymi na podstawie próby losowej.
1 2
Współczynnik ufności 1 - ą jest wielkością, którą można interpretować w następujący spo-
sób: jest to prawdopodobieÅ„stwo, że rzeczywista wartość parametru ¸ w populacji znajduje
się w wyznaczonym przez nas przedziale ufności. Im większa wartość tego współczynnika,
tym szerszy przedział ufności, a więc mniejsza dokładność estymacji parametru. Im mniejsza
wartość 1 - ą, tym większa dokładność estymacji, ale jednocześnie tym większe prawdopodo-
bieństwo popełnienia błędu. Wybór odpowiedniego współczynnika jest więc kompromisem
pomiędzy dokładnością estymacji a ryzykiem błędu. W praktyce przyjmuje się zazwyczaj
wartoÅ›ci zależne od parametru à - czyli odchylenia standartowego.
Odchylenie standardowe  klasyczna miara zmienności, obok średniej arytmetycznej naj-
częściej stosowane pojęcie statystyczne. Intuicyjnie rzecz ujmując, odchylenie standardowe
mówi, jak szeroko wartości jakiejś wielkości są rozrzucone wokół jej średniej. Im mniejsza
wartość odchylenia tym obserwacje są bardziej skupione wokół średniej. Odchylenie standar-
dowe jest pierwiastkiem kwadratowym z wariancji. Pojęcie odchylenia zostało wprowadzone
przez pioniera statystyki, Karla Pearsona w 1894 roku. Odchylenie standardowe zmiennej lo-
sowej oznacza się tradycyjnie przez à (małe greckie sigma) i definiuje jako pierwiastek kwa-
dratowy wariancji. Jest ono dane wzorem:
2
à = E(( X - E( X )) )
gdzie: E(X) jest wartością oczekiwaną X.
Metody estymacji parametrów modeli parametrycznych pozwalają na jednoczesne
wyznaczenie macierzy kowariancji uzyskanych parametrów. Miarą zmienności zmiennej
losowej w statystyce jest wariancja definiowana jako, średnia arytmetyczna kwadratów
odchyleń poszczególnych wartości zmiennej od jej wartości oczekiwanej. Gdy rozważamy
kilka zmiennych losowych ich zależność liniową określa funkcja kowariancji:
Cxy = Cov(x, y) = E{(x - E{x})( y - E{ y})}
gdzie: C  kowariancja pomiędzy zmiennymi x i y,
xy
E  wartość oczekiwana
Mając macierz kowariancji na głównej przekątnej możemy znalez
ć wartości wariancji
kolejnych parametrów modelu, a z nich pierwiastkując policzyć odchylenia standardowe. Jak
już powiedziano przedziały ufności wyznacza się na podstawie znajomości odchylenia
standartowego. I tak wyróżnić możemy: metodÄ™ 1 à gdzie, dla rozkÅ‚adu normalnego, 68,3 %
populacji próby losowej znajduje siÄ™ w zaÅ‚ożonych przedziaÅ‚ach ufnoÅ›ci, metodÄ™ 2 à gdzie
mamy 95,4 % i wreszcie metodÄ™ 3 à gdzie mamy 99,7 % próby w przedziale ufnoÅ›ci.
Modele modalne
W ocenie poprawności estymacji modeli modalnych stosuje szereg kryteriów.
Najpopularniejszym z nich jest MAC (Modal Assurance Criterion). Wyznaczenie go polega
na wyliczaniu relacji między dwoma wektorami modalnymi wg następującego wzoru:
2
*
{Vrtest}{VrES}
'
MACr,r ' =
* T *
({Vrtest}{Vrtest} )({VrES} {VrES} )
' '
Przyjmuje on wartość od 0 do 1. Wysoka wartość MAC wskazuje, że dwa porównywane
wektory modalne są podobne. Ma to duże znaczenie przy wyborze najlepszego modelu do
dalszych analiz oraz przy dostrajaniu modeli numerycznych do wyników eksperymentu.
Zgodność modeli eksperymentalnych i numerycznych ma kapitalne znaczenie w procesie
wirtualnego prototypowania, gdzie odpowiednie dostrojenie modelu MES jest podstawowym
warunkiem umożliwiającym przeprowadzenie poprawnych symulacji. W przypadku, gdy
model nie będzie dostrojony  zgodny z rzeczywistością wyniki obliczeń przeprowadzanych z
zastosowaniem takiego modelu nie będą wiarygodne. Pochodnymi kryterium MAC są
MACCo i CoMAC.
Do wskazania, które stopnie swobody mają największy wpływ na wartość MAC służy
kryterium MACCo (Modal Assurance Criterion Contribution). Pokazuje ono, które stopnie
swobody należy usunąć z modelu (wektorów modalnych), aby wartość MAC wzrosła do
zadanego poziomu. Pozwala to na wskazanie miejsc o najwyższej wrażliwości na zakłócenia i
eliminację błędnych pomiarów.
Aby porównać, w jaki sposób poszczególne stopnie swobody są ze sobą skorelowane posłu-
żyć należy się wskaz
nikiem CoMAC (Coordinated Modal Assurance Criterion). Sprawdza on
podobieństwo wektorów własnych obu modeli dla wszystkich PDW w kolejnych kierunkach
pomiarowych. Rysunek poniżej przedstawia przykładowe wartości wskaz
nika CoMAC  war-
tości bliskie 1 oznaczają dobrze skorelowane stopnie swobody.
Modele nieparametryczne
Przejdz
my teraz do oceny modeli nieparametrycznych. Podstawowym modelem
nieparametrycznym jest widmowa funkcja przejścia układu dynamicznego. Funkcją pochodną
do niej jest funkcja koherencji zwyczajnej. Funkcja ta informuje nas jaka część odpowiedzi
dynamicznej układu (użytej do estymacji WFP) została wywołana przez mierzone
wymuszenie układu (również użyte do estymacji WFP). Innymi słowy daje ona pogląd na to
na ile odpowiedz
układu wywołana była danym wymuszeniem. Funkcja koherencji
2
zwyczajnej oznaczana jest przez ł I jest wyliczana z następującego wzoru:
pq
gdzie: GXF wzajemna gęstość widmowa mocy pomiędzy odpowiedzią mierzoną w punkcie
pq
p a wymuszeniem zadawanym w punkcie q.
Funkcja koherencji zwyczajnej jest określana w dziedzinie częstotliwości I przyjmuje
wartości od 0 do 1. Jest ona równa 1 dla danej częstotliwości gdy zmierzona gęstość
widmowa mocy jest wywołana całkowicie i jedynie przez zmierzoną gęstość widmową mocy
wymuszenia. Wartość koherencji mniejsza niż 1 wskazuje, że zmierzona odpowiedz
jest
większa od tej wywołanej danym wymuszeniem. Może to być spowodowane działaniem
innych z
ródeł, szumem pomiarowym lub nieliniowościami w układzie. Gdy wartość funkcji
koherencji jest niska do poprawnej estymacji WFP potrzebnych jest znacznie więcej
uśrednień. Gdy wartość koherencji jest równa zero, odpowiedz
układu jest całkowicie
wywołana przez z
ródła inne inż. Mierzone. Generalnie koherencja może być uznawana za
miarę udziału szumu w pomiarach. Na rysunku poniżej pokazano przykładowy przebieg
funkcji koherencji dla danej WFP.
l e ft : 3 : Y : l e ft : 3 : Y
1 0
5
0
- 5
- 1 0
- 1 5
- 2 0
- 2 5
- 3 0
0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0
H z
l e ft : 3 : Y : l e ft : 3 : Y
1
0 . 9
0 . 8
0 . 7
0 . 6
0 . 5
0 . 4
0 . 3
0 . 2
0 . 1
0
0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0
H z
Na powyższym rysunku widać wyraz
nie, że wartości koherencji dla danego pomiaru są
odpowiednie dla zakresu częstotliwości 35  256 Hz. Poniżej tego pasma odpowiedz
układu
nie była powodowana jedynie przez mierzone wymuszenie. Rysunek ten pokazuje też inne
ważne zjawisko, mianowicie spadek wartości koherencji w częstotliwościach rezonansowych
i antyrezonansowych układu.
Istnieją też kryteria, które pozwalają na porównywanie modeli nieparametrycznych (np.
WFP) uzyskanych w różnych eksperymentach i badaniach symulacyjnych w celu wyłonienia
najlepszego modelu, bÄ…dz
znalezienia zgodności pomiędzy modelem eksperymentalnym (rze-
czywistym) i numerycznym, zbudowanym np. przy pomocy MES. Z grupy tych kryteriów
wybrano RVAC (Response Vector Assurance Criterion) i FRAC (Frequency Response
Assurance Criterion).
Pierwsze kryterium rozpatruje wektory wartości amplitud WFP dla wszystkich stopni swobo-
dy w danej częstotliwości. Następnie liczony jest MAC między wektorem stworzonym z WFP
z pierwszego pomiaru i wektorem WFP pochodzÄ…cych z innego pomiaru lub wysyntezowa-
nych przy pomocy modelu ES. Operacja ta powtarzana jest dla całego zakresu częstotliwości i
z uzyskanych wartości MAC tworzony jest wykres w dziedzinie częstotliwości. Wartości
RVAC bliskie 1 wskazują na duże podobieństwo WFP w danej częstotliwości.
-
d B ( m * s 2 / N )
/
RVAC(É) = MAC({E (É)},{U (É)})
test fem
FRAC wyliczany jest ze wzoru:
2
*
ES test
{H(É ) }{H(É ) }
i i
j j
FRAC( j) =
* *
ES ES test test
ëÅ‚ öÅ‚ ëÅ‚ öÅ‚
{H(É ) }{H(É ) } {H(É ) }{H(É ) }
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚
i i i i
j j j j
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
gdzie: {H(É ) } jest wektorem opisujÄ…cym WFP dla stopnia swobody j dla zakresu czÄ™stotli-
i j
woÅ›ci É . Dobrze skorelowane WFP dajÄ… FRAC bliski 1.
i
Ponieważ model elementów skończonych nieco zawyża sztywność obiektu, dlatego następuje
przesunięcie w górę częstotliwości PDW. Efektem tego mogą być niskie wartości FRAC
mimo faktycznego podobieństwa. Należy więc dobrać taką wartość współczynnika przesunię-
cia, WFP aby częstotliwości drgań własnych dla obu modeli się pokryły. Aby znalez
ć najlep-
szą zgodność WFP, ustala się pewien zakres zmienności współczynnika przesunięcia i spraw-
dza, która wartość daje najlepsze wyniki.
Zadania do samodzielnego wykonania:
1. Zamodeluj analogicznie jak w przykładzie układ z rysunku poniżej i przeprowadz
jego
symulację wymuszając go szumem losowym. Przebiegi czasowe zarówno wejścia (siła f )
1
jak i wyjścia (przemieszczenia x ) zapisz i dodaj do nich szum losowy o wartości amplitu-
2
dy = 20 % amplitudy sygnału zaszumianego. Następnie dokonaj estymacji parametrów
modeli typu AR, OE i PEM i porównaj te modele przy pomocy kryterium AIC. Skomentuj
wyniki.
Dane: M = 10, M = 7, C = 0.7, C = 0.9, C = 0.7, K = 120, K = 200, K = 170
1 2 1 2 3 1 2 3
2. Dla modelu z zadania 1 dokonaj estymacji WFP. Następnie zmień sztywności K1 i K2
o około 10 % i ponownie dokonaj syntezy WFP. Porównaj otrzymane charakterystyki przy
pomocy kryteriów FRAC i RVAC.
3. Zastosuj przebiegi wejścia (siła f ) jak i wyjścia (przemieszczenia x ) z zadania 1 do
1 2
estymacji funkcji koherencji zwyczajnej (oba przebiegi bez zakłóceń). Następnie dodaj do
przebiegu wyjścia szum losowy o amplitudzie 30 % amplitudy sygnału zaszumianego i
ponownie policz funkcję koherencji zwyczajnej. Czy 30 % szum na wyjściu spowodował
adekwatny spadek wartości funkcji koherencji?
4. Wygeneruj dwa zestawy 7 siedmioelementowych wektorów losowych. Dla wygenero-
wanych wektorów wyznacz macierz wskaz
ników MAC. Następnie wez
tylko pierwszy ze-
staw i policz autoMAC (MAC pomiędzy dwoma identycznymi zestawami wektorów). Na-
stępnie do pierwszego zestawu wektorów dodaj szum losowy o wartości około 10 % am-
plitudy zakłócanych wektorów i jeszcze raz policz MAC pomiędzy zestawem wektorów
bez i z zakłóceniem.