Rozpoznawanie parametrów modelu

Mammoth

2009-10-16 23:06:53
2009-10-16 23:06:53
2 komentarze
3059 odsłon

1 Wstęp
2 Opis problemu
3 Algorytm
4 Implementacja
5 Działanie programu
     5.1 Wyliczanie stałej g
     5.2 Bramka logiczna AND
     5.3 Przykład z „optymalizacją”
6 Podsumowanie

Wstęp

W tym artykule chciałem przedstawić mały program do rozpoznawania parametrów jakiegoś modelu na podstawie tego jak reaguje on na zadane sygnały wejściowe. Z polskiego na nasze. Mamy jakiś obiekt, oddziaływujemy na niego jakimiś sygnałami wejściowymi i obserwujemy jak na nie reaguje. Następnie podajemy jakiś hipotetyczny model (zestaw opisujących go równań) i na podstawie zebranych danych obliczmy współczynniki przy poszczególnych składowych równań (tu wkracza mój program). Dzięki temu dostajemy pełny model opisujący dany obiekt (oczywiście z pewną dokładnością). Może brzmi to jak jakiś żargon więc podam jakiś przykład do czego może to się przydać – wyliczanie stałych fizycznych na podstawie zebranych danych doświadczalnych, określanie równania krzywej najlepszego dopasowania na wykresach oraz – jeżeli się trochę postaramy – identyfikacja obiektu lub nawet optymalizacja równania.

A propos rysowania wykresów. Mój program do realizacji zadania posługiwał się będzie metodą tzw. „najmniejszych kwadratów” (brzmi znajomo). Jest to dobra i dość dokładana metoda.

Wiem, że może nie za ciekawie zaczyna się ten artykuł, ale postaram się opisywać wszystko w miarę przystępnie i przejrzyście.
Zadaje sobie również sprawę z tego, że poniższy program można zrealizować w 3 linijkach używając Matlab’a, Octave lub Exel’a. Jednak jest to mój pierwszy program który coś liczy i którego można wykorzystać do rozwiązywania realnych problemów więc postanowiłem się nim podzielić ;-) .
Bardzo proszę o zgłaszanie wszelkich błędów (merytorycznych jak i tych w kodzie programu) w komentarzach. Proszę tez o to by nie zaczynały się one od „Ty debilu..” lub „Nieuku..”. Zdaję sobie sprawę, że mój program nie jest niczym szczególnym oraz że mogę się mylić w niektórych miejscach. Najwyżej będziemy potem poprawiać.

Opis problemu

Jeżeli ktoś ma powyżej uszu matematykę, wie o co chodzi bądź woli skupić się na implementacji proponuje przejść do następnego podpunktu.

Jak wspomniałem większość rzeczy można opisać na zasadzie modelu mającego jakieś wejście (u) oraz wyjście (y). Znając model matematyczny jesteśmy w stanie przewidzieć co będzie na wejściu lub na odwrót mając pewne dane doświadczalne możemy określić równania.

Większość modeli można sprowadzić do równania macierzowego postaci:

(1) $y=au$

Muszę tu jednocześnie nadmienić że chodzi tu o modele liniowe, mój program nie oblicza równań dynamicznych (opisywanych całkami, bądź równaniami różniczkowymi).

Równanie (1) to można rozpisać:

(2) $y=\alpha \varphi (u)$

$\alpha$ – macierz współczynników. Tego szukamy.
$\varphi$ – klasa modelu. Jest to macierz prawdopodobnie opisujących model równań rozbitych na cząstki (wg. znaku plus). Ponadto równania te są pozbawione współczynników. Wygląd tej macierzy pokazuje poniższy przykład - równanie (4).
$y$ – macierz wyjść.
$u$ – macierz wejść

Przykład:

Mamy obiekt z trzema wejściami oraz dwoma wyjściami. Model matematyczny opisujący go może wyglądać następująco:

(3) $ \left{ \begin{array}{ll} y_1=\alpha _{11} u_1^2 + \alpha _{12} u_2 + \alpha _{13} u_3 / u_2 + \alpha _{14} \cdot 0 \ y_2=\alpha _{21} u_1 u_3 + \alpha _{22} u_3 + \alpha _{23} u_1^3 + \alpha _{24} (u_2 + u_3) \end{array} \right. $

Zapisując ten model w postaci takiej jak (2) będzie on wyglądał następująco:

(4) $ \left( \begin{array}{c} y_1 \ y_2 \end{array} \right) = \left( \begin{array}{cccc} \alpha _{11}& \alpha _{12} & \alpha _{13} & \alpha _{14} \ \alpha _{21}& \alpha _{22} & \alpha _{23} & \alpha _{24} \end{array} \right) \cdot \left( \begin{array}{cc} u_1^2 & u_1u_3 \ u_2 & u_3 \ u_3/u_2 & u_1^3 \ 0 & u_2+u_3 \end{array} \right) $

Naszym głównym problemem jest wyliczenie α. Nie możemy po prostu podać, że jest to

(5) $\alpha = \varphi ^{-1} (u) \cdot y$

bo musimy uwzględnić ilość pomiarów (im więcej tym dokładniej określimy model). Stąd też wynika jeden z najważniejszych problemów. Nasz model będzie zawsze tylko pewnym przybliżeniem rzeczywistości. Spowodowane jest to m.in. przez błędy pomiarowe, ale gł. przez to, że jest to tylko matematyczny model uwzględniający tylko podstawowe parametry. Dlatego nasz model nie będzie z reguły w stanie wyliczyć tych samych wartości wyjściowych (y) zebranych podczas doświadczenia przy danych wartościach wejściowych (u). Będzie on natomiast wyliczał pewne mniej lub bardziej rozbieżne wyjścia.

(6) $\bar{y} = \alpha \varphi (u)$

Różnica między y a y' jest podstawą wielu tzw. kryteriów najlepszego dopasowania Q(α). Jak już wspomniałem wcześniej mój program realizuje tylko jeden (najmniejszych kwadratów). Poniżej pokarzę jak można dzięki niemu wyliczyć równanie na α* (α najlepiej dopasowane). Na tej podstawie będzie można wyliczyć inne pod warunkiem, że będą one liniowe – będą miały pochodną (egh…) w każdym punkcie.

Ogólnie kryterium dopasowania wygląd tak:

(7) $Q(\alpha *) = min Q(\alpha )$

By było ono minimalne narzucamy warunki:

(8) $\frac{ \partial Q( \alpha ) }{ \partial \alpha } = 0$
(9) $\frac{ \partial Q( \alpha )}{ \partial \alpha ^2 } > 0$

Tzn. pierwsza pochodna równa zero oraz druga większa od zera. Czyli określamy minimum lokale funkcji. Kryterium najmniejszych kwadratów wygląda następująco (rozpisane):

(10) $Q( \alpha ) = \displaystyle\sum_{i=1}N {(y_i - \bar{y}i )2}= \sum{i=1}N {[y_i - \alpha \varphi(u_i)]2}= \sum_{i=1}N {y_i2} - 2\alpha \sum_{i=1}N {y_i \varphi (u_i)} + \alpha 2 \sum_{i=1}N {\varphi 2 (u_i)}$

Gdzie N to ilość pomiarów. Mówiąc opisowo metoda najmniejszych kwadratów polega na tym by suma kwadratów poszczególnych odchyleń od wartości rzeczywistej była jak najmniejsza.

Stosując warunek (8) otrzymujemy:

(11) $ \frac {\partia [ \sum_{i=1}N {y_i2}]}{\partial \alpha } - \frac {\partia [2\alpha \sum_{i=1}N {y_i \varphi (u_i)}]} {\partial \alpha } + \frac {\partia [\alpha 2 \sum_{i=1}N {\varphi 2 (u_i)}]} {\partial \alpha } = \ 0 - 2 \sum_{i=1}N {y_i \varphi (u_i)} + 2 \alpha \sum_{i=1}N { \varphi ^2(u_i)} = \ \alpha \sum_{i=1}N { \varphi 2(u_i)} - \sum_{i=1}^N {y_i \varphi (u_i)} = 0$

Stąd:

(12) $ \alpha * = \frac { \sum_{i=1}N {y_i \varphi (u_i)} } { \sum_{i=1}N { \varphi ^2(u_i)} }$

Stosując warunek (9) otrzymamy:

(13) $ \frac{ \partial Q( \alpha )}{ \partial \alpha 2 } > 0 \Rightarrow \sum_{i=1}N {\varphi ^2 (u_i)} > 0$

Czyli przynajmniej jeden $\varphi(u) \not= 0$ . Musimy o to zadbać przy korzystaniu z programu.

W ten sposób otrzymaliśmy równanie na α (przynajmniej najlepiej dopasowane pod względem kryterium „najmniejszych kwadratów”).

Algorytm

Zamieniając równanie (12) na postać macierzową otrzymamy:

(14) $ \alpha * = [ \sum_{i=1}N {\varphi (u_i) \cdot \varphi T (u_i)} ]{-1} \cdot \sum_{i=1}N {y_i \varphi (u_i)}$

$\alpha *$ – macierz szukanych współczynników
$\varphi$ – klasa modelu (macierz równań)
$u$ – wektor danych wejściowych
$y$ – wektor danych wyjściowych

Od razu widzimy jak nasz algorytm powinien wyglądać oraz jakich danych potrzebuje.

Słowna postać algorytmu:

Potrzebne dane:
Macierz równań Fi,
Macierz danych wejściowych u,
Macierz danych wyjściowych y.
Potrzebne zmienne:
Macierz liczbowa f,
Macierz liczbowa fT,
Macierz liczbowa b (reprezentująca drugą sigmę w równaniu (14)),
Macierz liczbowa FI (reprezentująca pierwszą sigmę w równaniu (14)).

Dla każdego z pomiarów (od i=1 do N)
a) wylicz fi(ui)
b) wypełnij fiT(ui)
c) dosumuj do macierzy b wartość fi(ui)*yi
d) dosumuj do macierzy FI wartość fi(ui)*fiT(ui)
obróć macierz FI
pomnóż FI przez b
zwróć wynik w postaci FI

Implementacja

Podczas implementacji pojawiły się dwa zasadnicze problemy:

wygodna reprezentacja typu danych Macierz zdolnego przechować liczby jak i równania,
wygodna reprezentacja typu danych Równanie.
Do rozwiązania tych problemów skorzystałem z udogodnień języka C++ takich jak przeciążanie operatorów oraz szablony.

W pierwszym przypadku stworzyłem szablon macierzy. Przyznam, że zwyczajnie przerobiłem jeden z projektów zrobionych na laboratorium z programowania ;-). Projekt ten zawierał już mnożenie, dodawanie, odwracanie, potęgowanie, … Wystarczyło zamienić typ double na TYP i tyle.

W drugim przypadku było nieco gorzej. Z równania (14) widać, że macierz równań fi(u) jest wymagana w trzech miejscach algorytmu. Oczywistym jest, że dla każdej z iteracji musi być ona wyliczona oddzielnie. Postanowiłem więc, że Fi będzie macierzą przechowującą równania RPN. W każdej iteracji zostanie obliczony jej odpowiednik liczbowy (w zależności od ui). Tak stworzona macierz zostanie użyta do dalszych obliczeń.

Sam algorytm przedstawia się następująco:

Matrix<double> Model::Identify()
{
	unsigned m = _eq.Rows(), n = _y.Rows();

	Matrix<double> fI(m, m), b(m, n),
	               f(m, n), ft(n, m);

	fI.O(); b.O();  //wyzerowanie macierzy

	for (unsigned i=0; i<_u.Rows(); i++)
	{
		//oblicza fi(ui)
		for (unsigned x=0; x<m; x++)
			for (unsigned y=0; y<n; y++)
				f[x][y] = _eq[x][y].Solve(_u[i]);
		//wypelnia fi^T(ui)
		for (unsigned x=0; x<m; x++)
			for (unsigned y=0; y<n; y++)
				ft[y][x] = f[x][y];
		//oblicza yi*fi(ui)
		for (unsigned x=0; x<m; x++)
			for (unsigned y=0; y<n; y++)
				b[x][y] += (_y[y][i]*f[x][y]);
		//
		fI += (f*ft); //Wiem, że nie jest to najlepsze
                              //rozwiązanie bo powoduje narzuty
                              //pamięciowe, ale przynajmniej dobrze
                              //(czytelnie) wygląda ;-)
	}

	//oblicza fi(u)*fi^T(u)*y*f(u)
	//fI == Sigma(fi(ui)*fi^T(ui))
	//b == Sigma(yi*f(ui))
	(fI^=-1)*=b;

	return fI;
}

Cały projekt zorganizowany jest w następujący sposób:

Uznałem, że nie ma sensu publikowania całego kodu w artykule, można go jednak pobrać RozpoznawanieParametrowModelu.zip.

Działanie programu

Program działa w trybie tekstowym, a jego obsługa jest bardzo łatwa. Poniżej demonstruje kilka listingów z jego działania.

Wyliczanie stałej g

Doświadczenie polegało na zrzuceniu obiektu z pewnej wysokości. W kolejnych sekundach mierzyliśmy przebytą drogę.

$u_1 - t[s]$	1	2	3	4	5	6
$y_1 - s[m]$	8.49	20.05	50.65	72.19	129.85	179.56

Oczywiście korzystamy ze wzoru:

$s = \frac {1} {2} g t ^2$

Którego może przekształcić dla naszych potrzeb:

$y_1 = \alpha _1 u_12 = \alpha _1 \cdot \frac {1} {2} \cdot u_12$

Zastosujemy więc model z parametrem 0.5. Pozwoli to na bezpośrednie wyliczenie szukanej stałej.

Ilosc pomiarow: 6
Ilosc wejsc: 1
u1: 1 2 3 4 5 6
Ilosc wyjsc: 1
y1: 8.49 20.05 50.65 72.19 129.85 179.56

u
| 1 |
| 2 |
| 3 |
| 4 |
| 5 |
| 6 |

y
| 8.49 20.05 50.65 72.19 129.85 179.56 |

Liczba modeli do zbadania: 1

 ===> Model #1 ===>
Dl. najdluzszego rownania: 1
y1 = 0.5*u1^2

| {0.5 u1 2 ^ * } |

Dopasowanie:
| 10.0308 |

y':
5.01538 e=3.47462
20.0615 e=-0.0115209
45.1384 e=5.51158
80.2461 e=-8.05608
125.385 e=4.46549
180.554 e=-0.993688

I tak oto otrzymaliśmy krzywą najlepszego dopasowania.

Wyjaśnienia może wymagać pole „Dl. najdluzszego rownania”. Tą wartością podpowiadamy programowi jak szeroka powinna być macierz równań (ile składowych względem plusa posiada najdłuższe równanie). Jeżeli jakieś równanie jest krótsze od podanej wartości wtedy pozostałe komórki macierzy inicjowana są zerami.

Bramka logiczna AND

Zwyczajne określenie wzoru opisującego bramkę AND. Należy tutaj pamiętać, że jeżeli spodziewamy się jakiejś interakcji między jakimiś wejściami możemy to uwzględnić w równaniu dodając iloczyn owych wejść. W naszym przypadku wpisaliśmy na oślep równanie

$y_1 = u_1 + u_2 + u_1u_2$

Ilosc pomiarow: 4
Ilosc wejsc: 2
u1: 1 1 0 0
u2: 1 0 1 0
Ilosc wyjsc: 1
y1: 1 0 0 0

u
| 1 1 |
| 1 0 |
| 0 1 |
| 0 0 |

y
| 1 0 0 0 |

Liczba modeli do zbadania: 1

 ===> Model #1 ===>
Dl. najdluzszego rownania: 3
y1 = u1+u2+u1*u2

| {u1 } |
| {u2 } |
| {u1 u2 * } |

Dopasowanie:
| 1.11022e-016 |
| 1.11022e-016 |
| 1 |

y':
1 e=2.22045e-016
1.11022e-016 e=-1.11022e-016
1.11022e-016 e=-1.11022e-016
0 e=0

W efekcie dostaliśmy dopasowanie $y_1 = 0 \cdot u_1 + 0 \cdot u_2 + 1 \cdot u_1u_2 = u_1u_2$ , czyli „zidentyfikowaliśmy” obiekt.

Należy tutaj dodać że liczba postaci 1.11022e-016 jest niezmiernie mała i bliska zeru. Jest to wynik ograniczeń reprezentacji zmiennoprzecinkowej w komputerze.

Przykład z „optymalizacją”

W zwyczajnym Exel’u wklepałem kilka danych wejściowych i podając pewną formułę wygenerowałem jakieś dane wyjściowe. Powiedzmy, że jest to rzeczywiste równanie opisujące obiekt.

$\left\{ \begin{array}{ll} y_1= u_1 + u_2 + u_3 \\ y_2= u_1 u_3 \end{array} \right.$

$u_1$	1	2	3	1	3
$u_2$	1	2	3	2	2
$u_3$	1	2	3	3	1
$y_1$	3	6	9	6	6
$y_2$	1	4	9	3	3

Ilosc pomiarow: 6
Ilosc wejsc: 3
u1: 0 1 2 3 1 3
u2: 0 1 2 3 2 2
u3: 0 1 2 3 3 1
Ilosc wyjsc: 2
y1: 0 3 6 9 6 6
y2: 0 1 4 9 3 3

u
| 0 0 0 |
| 1 1 1 |
| 2 2 2 |
| 3 3 3 |
| 1 2 3 |
| 3 2 1 |

y
| 0 3 6 9 6 6 |
| 0 1 4 9 3 3 |

Liczba modeli do zbadania: 1

 ===> Model #1 ===>
Dl. najdluzszego rownania: 3
y1 = u1+u2+u3
y2 = u1*u3

| {u1 } {u1 u3 * } |
| {u2 } {0 } |
| {u3 } {0 } |

Dopasowanie:
| 6.66134e-016 1 |
| 3 -2 |
| 7.10543e-015 1 |

y':
0 e=0 0 e=0
3 e=-8.88178e-016 1 e=0
6 e=-1.77636e-015 4 e=0
9 e=-1.77636e-015 9 e=0
6 e=-7.99361e-015 3 e=0
6 e=5.32907e-015 3 e=0

Czyli w zasadzie dopasowanie wyszło nam:

$\left\{ \begin{array}{ll} y_1= 3u_2 \\ y_2= u_1 u_3 \end{array} \right.$

W efekcie mam optymalizację równania, jednak tylko dla tak specyficznego zestawu danych wejściowych i wyjściowych.

Podsumowanie

Trochę tandetnie to wyszło – tyle wzorów i pitolenia, ale cóż ciężko było to jakoś inaczej przedstawić. Mam nadzieję, że program się spodoba – a co ważniejsze – sprowokuje pozostałych do pisania innych tego typu numerycznych dziwadeł. Wiem, że program ma pewne ograniczenia i niedociągnięcia - więc nie znęcajcie się nade mną w komentarzach ;-). Wniosek z powyższego artykułu jest taki, że rozwiązując nawet najbardziej zwyczajny problem (rysowanie wykresu) można się sporo nauczyć i ciekawego dowiedzieć. Co więcej, znienawidzona przez wszystkich – i mnie też – matma może się do czegoś przydać ;-)

Załączniki

kod źródłowy RozpoznawanieParametrowModelu.zip.

Autor: Mammoth (maf.michal@gmail.com)

Rozpoznawanie parametrów modelu

Mammoth

Wstęp

Opis problemu

Algorytm

Implementacja

Działanie programu

Wyliczanie stałej g

Bramka logiczna AND

Przykład z „optymalizacją”

Podsumowanie

2 komentarzy

Praca dla programistów

Forum dyskusyjne

Sprawy administracyjne

O nas

Skontaktuj się z nami