MATLAB - Aproksymacja algorytmem LS

Kontynuuje omawianie metody LS do aproksymacji danych pomiarowych. Tym razem na warsztat idzie algorytm LS w postaci macierzowej, który umożliwia zastosowaniem modeli wyższych rzędów. Zobaczmy jak dopasować parabolę do danych układających się w ten właśnie kształt.

W poprzednim wpisie omówiłem estymator Least Squares (LS) w ujęciu algebraicznym – wyprowadziłem wzory na współczynniki prostej, aproksymującej dane pomiarowe. Jest to jednak przypadek szczególny, kiedy mamy do czynienia z liniową relacją między argumentem u, a wartościami y. Gdyby jednak relacja u-y była nieliniowa, model w postaci prostej byłby błędnym wyborem – należałoby zastosować model w postaci krzywej (wielomianu) wyższego rzędu. W takim przypadku nie da się w prosty sposób wyprowadzić analitycznie wzorów np. na równanie paraboli. Dużo prościej, i tak robi się w praktyce, jest zastosować estymator LS w postaci macierzowej.

Algorytm / estymator LS

Algorytm taki jest również zastosowany w pakiecie funkcji z toolboxu curve fitting MATLABa. Jakkolwiek można stosować te funkcje w sposób automatyczny, to dobrze jest również wiedzieć jak one działają i jakie cechy ma zastosowany tam estymator. Dzisiaj napiszemy własny program realizujący aproksymację danych pomiarowych wielomianem dowolnego stopnia.

Na potrzeby opisu estymatora LS w postaci macierzowej, przyjmijmy, że mamy do czynienia z obiektem, na wejściu którego mamy sygnał wejściowy ui, a na wyjściu yi (patrz analogia do poprzedniego wpisu).

Narysujmy relacje pomiędzy sygnałem wyjściowym, a wejściowym. W tym przypadku jest ona nieliniowa:

Na potrzeby dalszych obliczeń, zapiszmy wszystkie wartości/próbki sygnału wyjściowego w postaci wektora:

Ten wektor będzie nam potrzebny pod koniec tych rozważań, gdyż to właśnie z nim będziemy porównywać nasz model, chcąc uzyskać jak największą zbieżność pomiędzy wynikami rzeczywistymi (z obiektu) a z modelu.

Z powyższego rysunku widać wyraźnie, że zastosowanie modelu w postaci równania prostej byłoby w tym przypadku błędem. Załóżmy więc, że ogólnie chcemy dopasować model w postaci wielomianu K-tego rzędu:

gdzie:

- i-ta próbka sygnału wyjścia (odpowiedzi) modelu, i=1,…,N.

– i-ta próbka sygnału wejściowego (wymuszenia) obiektu i=1,…,N.

- poszukiwane współczynniki modelu, j=0,…,K.

Ponieważ zarejestrowane sygnały yi i ui mają po N próbek, to powyższe równanie można zapisać jako układ równań, gdzie każde z nich odpowiada kolejnej próbce wymuszenia:

Pogrupujmy zmienne w postaci macierzy:

Możemy teraz zapisać powyższe równania w zwartej, macierzowej postaci:

Gdzie ym to odpowiedz dopasowywanego modelu, do danych pomiarowych zawartych w macierzy U (i wektorze y), za pomocą estymowanych (metodą LS) współczynników b. W jaki sposób obliczyć współczynniki naszego modelu? Podobnie jak miało to miejsce w metodzie algebraicznej (patrz poprzedni wpis), należy utworzyć kryterium Q, będące liczbową miarą „odległości” pomiędzy danymi pomiarowymi, a dopasowywanym modelem, i znaleźć jego minimum ze względu na współczynniki b:

Nie będę tutaj zamieszczał obliczania pochodnej w przypadku takiego równania macierzowego – takie wyprowadzenie możecie z znaleźć w literaturze. Końcowy wzór na poszukiwane współczynniki modelu wynosi:

Otrzymaliśmy uniwersalny, prosty we formie i łatwy do zastosowania wzór na współczynniki modelu dowolnego rzędu, który aproksymuje dane pomiarowe.

Algorytm Least Squares - rząd modelu

Pozostaje jeszcze odpowiedzieć na pytanie, jak ustalić rząd modelu, który chcemy dopasować do danych? Bardzo prosto, przez odpowiednie ukształtowanie macierzy U. Jeżeli chcemy otrzymać model pierwszego rzędu (prosta) pozostawiamy w tej macierzy tylko te kolumny, które odpowiadają argumentowi u i kolumnę jedynek (wyraz wolny):

Jeżeli natomiast chcemy dopasować parabolę czyli wielomian drugiego rzędu, to macierz U:

Analogicznie postępujemy w przypadku modeli wyższych rzędów. I to wszystko! Stosując jeden wzór oraz odpowiednio kształtując macierz U. możemy w łatwy sposób obliczyć współczynniki modelu aproksymującego dane pomiarowe. Najlepiej zobaczyć jak to działa na przykładzie (dane pomiarowe można pobrać TUTAJ).

clear all; close all; clc

load cz_pot2_nlin
N=length(u);

% wykres danych pomiarowych:
figure;
plot(u,y,'x','LineWidth',2);

%% Estymator LS w postaci algebraicznej (model liniowy)

b1=(N*sum(u.*y)-sum(u)*sum(y))/(N*sum(u.^2)-sum(u).^2)
b0=(sum(y)-b1*sum(u))/N

% przyjęty model:
ym=b1*u+b0;

% rysunek
hold on; plot(u, ym,'r','LineWidth',2); 
xlabel('u'); ylabel('y'); grid on;


%% Estymator LS w postaci macierzowej (model liniowy)

% macierz wejść:
U = [u ones(size(u))];

% estymator LS w postaci macierzowej:
b1 = inv(U'*U) * U' * y;

% model liniowy:
ym1 = U * b1;

% rysunki 
plot(u,ym1,'g','LineWidth',2);


%% Estymator LS w postaci macierzowej (model drugiego rzędu)

% macierz wejść:
U = [u.^2 u ones(size(u))];

% estymator LS w postaci macierzowej:
b2 = inv(U'*U) * U' * y;

% model nieliniowy:
ym2 = U * b2;

% rysunki 
plot(u,ym2,'b','LineWidth',2);
legend('dane pomiarowe', 'model liniowy (algebraicznie)', 'model liniowy (macierzowo)', 'model 2go rzędu (macierzowo)', 'location','northwest')
xlabel('u'); ylabel('y');

%% Porównanie residuł

figure;
bar(u, y-ym1);
hold on; bar(u, y-ym2);
xlabel('u'); ylabel('y-ym'); grid on
legend('model liniowy', 'model 2go rzędu')

Z obydwu rysunków, a w szczególności analizy wartości residuów wynika, że w omawianym przypadku, model drugiego rzędu jest zdecydowanie lepszym wyborem do aproksymacji danych pomiarowych.

Teraz jesteście gotowi do modelowania o dowolnym rzędzie!

(Visited 562 times, 1 visits today)

Dodaj komentarz

Twój adres e-mail nie zostanie opublikowany. Wymagane pola są oznaczone *