Akademia Górniczo-Hutnicza
Katedra Robotyki i Mechatroniki
Identyfikacja i analiza sygnałów
Laboratorium 11
Ocena jakości estymatorów
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 wskaź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ść:
)
ln(
2
2
L
k
AIC
−
−
=
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
PRZYKŁ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
1
, X
2
, ..., X
n
). Przedziałem ufności (θ - θ
1
, θ + θ
2
) o współczynniku uf-
ności 1 - α nazywamy taki przedział (θ - θ
1
, θ + θ
2
), który spełnia warunek:
P(θ
1
< θ < θ
2
) = 1 − α
gdzie: θ
1
i θ
2
są funkcjami wyznaczonymi na podstawie próby losowej.
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
X
E
X
E
−
=
σ
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:
{ }
(
)
{ }
(
)
{
}
y
E
y
x
E
x
E
y
x
Cov
C
xy
−
−
=
=
)
,
(
gdzie: C
xy
– kowariancja pomiędzy zmiennymi x i y,
E – wartość oczekiwana
Mając macierz kowariancji na głównej przekątnej możemy znaleźć 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
*
'
'
,
ES
r
T
ES
r
test
r
test
r
ES
r
test
r
r
r
V
V
V
V
V
V
MAC
=
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ę wskaź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 wskaźnika CoMAC – war-
tości bliskie 1 oznaczają dobrze skorelowane stopnie swobody.
Modele nieparametryczne
Przejdź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 odpowiedź układu wywołana była danym wymuszeniem. Funkcja koherencji
zwyczajnej oznaczana jest przez
γ
pq
2
I jest wyliczana z następującego wzoru:
gdzie: GXF
pq
wzajemna gęstość widmowa mocy pomiędzy odpowiedzią mierzoną w punkcie
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 odpowiedź jest
większa od tej wywołanej danym wymuszeniem. Może to być spowodowane działaniem
innych ź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, odpowiedź układu jest całkowicie
wywołana przez ź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.
0
5 0
1 0 0
1 5 0
2 0 0
2 5 0
3 0 0
- 3 0
- 2 5
- 2 0
- 1 5
- 1 0
- 5
0
5
1 0
H z
d
B
(m
*s
-
2
/N
)
l e f t : 3 : Y : l e f t : 3 : Y
0
5 0
1 0 0
1 5 0
2 0 0
2 5 0
3 0 0
0
0 . 1
0 . 2
0 . 3
0 . 4
0 . 5
0 . 6
0 . 7
0 . 8
0 . 9
1
H z
/
l e f t : 3 : Y : l e f t : 3 : Y
Na powyższym rysunku widać wyraźnie, że wartości koherencji dla danego pomiaru są
odpowiednie dla zakresu częstotliwości 35 – 256 Hz. Poniżej tego pasma odpowiedź 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ądź 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.
RVAC(
ω
) = MAC({E
test
(
ω
)},{U
fem
(
ω
)})
FRAC wyliczany jest ze wzoru:
( )
{
}
( )
{
}
( )
{
}
( )
{
}
( )
{
}
( )
{
}
=
*
*
2
*
)
(
test
j
i
test
j
i
ES
j
i
ES
j
i
test
j
i
ES
j
i
j
H
H
H
H
H
H
FRAC
ω
ω
ω
ω
ω
ω
gdzie: {H(
ω
i
)
j
} jest wektorem opisującym WFP dla stopnia swobody j dla zakresu częstotli-
wości
ω
i
. Dobrze skorelowane WFP dają FRAC bliski 1.
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 znaleźć 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 przeprowadź jego
symulację wymuszając go szumem losowym. Przebiegi czasowe zarówno wejścia (siła f
1
)
jak i wyjścia (przemieszczenia x
2
) zapisz i dodaj do nich szum losowy o wartości amplitu-
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
1
= 10, M
2
= 7, C
1
= 0.7, C
2
= 0.9, C
3
= 0.7, K
1
= 120, K
2
= 200, K
3
= 170
2. Dla modelu z zadania 1 dokonaj estymacji WFP. Następnie zmień sztywności K
1
i K
2
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
1
) jak i wyjścia (przemieszczenia x
2
) z zadania 1 do
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 wskaźników MAC. Następnie weź 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.